text:
function [] = exportInputsMustToGAMS(model, minFluxesW, maxFluxesW, constrOpt, inputFolder)
% This function export all the inputs needed to run the GAMS functions to
% find first order Must Sets (`MustU`, `MustL`). The inputs will
% be stored in inputFolder. Some inputs will be exported using GDXMRW and
% others will be exported as simple .txt files.
%
% USAGE:
%
%         exportInputsMustToGAMS(model, minFluxesW, maxFluxesW, constrOpt, inputFolder)
%
% INPUTS:
%    model:             (structure) a metabolic model with at least the
%                       following fields:
%
%                         * .rxns - Reaction IDs in the model
%                         * .mets - Metabolite IDs in the model
%                         * .S -    Stoichiometric matrix (sparse)
%                         * .b -    RHS of `Sv = b` (usually zeros)
%                         * .c -    Objective coefficients
%                         * .lb -   Lower bounds for fluxes
%                         * .ub -   Upper bounds for fluxes
%    minFluxesW:        (double array of size `n_rxns x 1`) minimum fluxes
%                       for each reaction in the model for wild-type
%                       strain. This can be obtained by running the
%                       function `FVAOptForce`. E.g.: `minFluxesW = [-90; -56];`
%    maxFluxesW:        (double array of size `n_rxns x 1`) maximum fluxes
%                       for each reaction in the model for wild-type
%                       strain. This can be obtained by running the
%                       function `FVAOptForce`. E.g.: `maxFluxesW = [90; 56];`
%    constrOpt:         (structure) structure containing additional
%                       contraints. Include here only reactions whose flux
%                       is fixed, i.e., reactions whose lower and upper
%                       bounds have the same value. Do not include here
%                       reactions whose lower and upper bounds have
%                       different values. Such contraints should be defined
%                       in the lower and upper bounds of the model. The
%                       structure has the following fields:
%
%                         * .rxnList - Reaction list (cell array)
%                         * .values -  Values for constrained reactions
%                           (double array). E.g.: `struct('rxnList', {{'EX_gluc', 'R75', 'EX_suc'}}, 'values', [-100, 0, 155.5]');`
%    inputFolder:       (string) Folder where inputs for GAMS function
%                       will be stored
%
% OUTPUTS:
%    Reactions.txt:     (file) File containing the identifiers for
%                       reactions
%    Metabolites.txt:   (file) File containing the identifiers for
%                       metabolites
%    Constrains.txt:    (file) File containing the identifiers for
%                       constrained reactions
%    MtoG.gdx:          (file) File containing the parameters which
%                       will be read by GAMS (lower bounds, upper bounds,
%                       stoichiometrix matrix `S`, minimum and maximun fluxes
%                       for each reaction in the previous step of FVA, and
%                       the values for contrained reactions)
%
% .. Author: - Sebastian Mendoza, May 30th 2017, Center for Mathematical Modeling, University of Chile, snmendoz@uc.cl

if nargin < 6 %input handling
    error('OptForce: All inputs must be specified when running exportInputsMustToGAMS')
end

%Create a temporaty folder for inputs
if ~exist(inputFolder, 'dir')
   mkdir(inputFolder);
end
current = pwd;
cd(inputFolder);

%Export Sets
exportSetToGAMS(model.rxns, 'Reactions.txt')
exportSetToGAMS(model.mets, 'Metabolites.txt')
exportSetToGAMS(constrOpt.rxnList, 'Constraints.txt')

%Export parameters
s.name = 's';
s.val = full(model.S);
s.type = 'parameter';
s.form = 'full';
s.uels = {model.mets',model.rxns'};

basemin.name = 'basemin';
basemin.val = minFluxesW;
basemin.type = 'parameter';
basemin.form = 'full';
basemin.uels = model.rxns';

basemax.name = 'basemax';
basemax.val = maxFluxesW;
basemax.type = 'parameter';
basemax.form = 'full';
basemax.uels = model.rxns';

lb.name = 'lb';
lb.val = model.lb';
lb.type = 'parameter';
lb.form = 'full';
lb.uels = model.rxns';

ub.name = 'ub';
ub.val = model.ub';
ub.type = 'parameter';
ub.form = 'full';
ub.uels = model.rxns';

b.name = 'b';
b.val = constrOpt.values;
b.type = 'parameter';
b.form = 'full';
b.uels = {'EX_gluc','R75','EX_suc'};

%Using GDXMRW to export inputs
fileName = ['MtoG' setType];
wgdx(fileName, s, basemin, basemax, lb, ub, b)

cd(current);

end

author: opencobra

repo: cobratoolbox


text:
function dicomrt_isosurface(inputdose,dose_xmesh,dose_ymesh,dose_zmesh,levels,norm,VOI,voi2plot)
% dicomrt_isosurface(inputdose,dose_xmesh,dose_ymesh,dose_zmesh,levels,norm,VOI,voi2plot)
%
% Show isosurface for 3D dose distributions.
%
% inputdose contains the 3D dose distribution
% inputdose can be an rtplan generated or a monte carlo generated 3D dose distribution or a dose difference.
%
% dose_xmesh,dose_ymesh,dose_zmesh are x-y-z coordinates of the center of the dose-pixel
% levels is a vector containing the values of the % dose difference to be displayed in 3D
% norm is the dose normalization level:
%   1.   norm =0  no normalization is carried out
%   2.       =-1 same as (1) but matrix is treated as dose difference (transparency display)
%   3.       >0  doses are normalized to norm (100%)
% VOI is an OPTIONAL cell array which contain the patients VOIs as read by dicomrt_loadVOI
% voi2plot is an OPTIONAL vector pointing to the number of VOIs to be displayed
%
% Example:
%
% dicomrt_isosurface(A,xmesh,ymesh,zmesh,[100 90],60,A_voi,[1 2 3])
%
% overlay 2 isosurfaces for matrix A at the 100% and 90% levels. Dose is normalised to 60Gy (100%).
% The figure also shows VOIs # 1 2 and 3.
%
% See also dicomrt_loaddose, dicomrt_loadmcdose, dicomrt_contourslice, dicomrt_dosediff
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

% Check number of argument and set-up some parameters and variables
error(nargchk(6,8,nargin))

% Check case and set-up some parameters and variables
[dose_temp,type_dose,doselabel,PatientPosition]=dicomrt_checkinput(inputdose);
dose=dicomrt_varfilter(dose_temp);

% Check normalization
if norm>0
    dose=dose./norm.*100;
end

% Initialize plot parameters:
color = char('r','g','b','c','m','y','k','w');
line = char('-','--',':','-.');
marker = char('+','o','*','x','s','d','^');
width = [0.5,1.0,1.5,2.0,2.5,3.0];
alphal=[0.1:0.1:1.0];

% Set isolevels display
[contourcolor,doseareas,transl]=dicomrt_setisolevels(norm,dose_temp);

% Fetching legend information
legend_matrix=char(num2str(levels'));

% Rendering 3D dose
handle=figure;
set(handle,'Name',['dicomrt_isosurface: ',inputname(1)]);
grid on;
hold

if exist('VOI')~=1 % No VOI is defined
    if length(levels)==1
        [n,p] = histc(levels,doseareas);
        handle = patch(isosurface(dose_xmesh,dose_ymesh,dose_zmesh,dose,levels));
        isonormals(dose_xmesh,dose_ymesh,dose_zmesh,dose,handle);
        set(handle,'FaceColor',contourcolor(p),'EdgeColor','None','FaceAlpha',0.5);
        view(30,30);
        legend(handle,[int2str(levels)]);
        xlabel('X axis (cm)','Fontsize',12)
        ylabel('Y axis (cm)','Fontsize',12)
        zlabel('Z axis (cm)','Fontsize',12)
    elseif length(levels)~=1
        for i=1:length(levels)
            [n,p] = histc(levels(i),doseareas);
            if p==0
                p=1;
            end
            handle = patch(isosurface(dose_xmesh,dose_ymesh,dose_zmesh,dose,levels(i)));
            isonormals(dose_xmesh,dose_ymesh,dose_zmesh,dose,handle);
            set(handle,'FaceColor',contourcolor(p),'EdgeColor','None','FaceAlpha',transl(p));
        end
        legend(legend_matrix);
        view(30,30);
        xlabel('X axis (cm)','Fontsize',12)
        ylabel('Y axis (cm)','Fontsize',12)
        zlabel('Z axis (cm)','Fontsize',12)
    end
    title(['Isosurface plot for: ',inputname(1)],'interpreter','none','FontSize',16);
elseif exist('VOI')==1 & exist('voi2plot')==1 % VOI was defined
    % Set VOI display parameters
    colorselect =7;
    lineselect =1;
    linewidth =1;
    if length(levels)==1
        [n,p] = histc(levels,doseareas);
        if p==0
            p=1;
        end
        handle = patch(isosurface(dose_xmesh,dose_ymesh,dose_zmesh,dose,levels));
        isonormals(dose_xmesh,dose_ymesh,dose_zmesh,dose,handle);
        set(handle,'FaceColor',contourcolor(p),'EdgeColor','None','FaceAlpha',0.5);
        view(30,30);
        legend(handle,[int2str(levels)]);
        xlabel('X axis (cm)','Fontsize',12)
        ylabel('Y axis (cm)','Fontsize',12)
        zlabel('Z axis (cm)','Fontsize',12)
    elseif length(levels)~=1
        for i=1:length(levels)
            [n,p] = histc(levels(i),doseareas);
            if p==0
                p=1;
            end
            handle = patch(isosurface(dose_xmesh,dose_ymesh,dose_zmesh,dose,levels(i)));
            isonormals(dose_xmesh,dose_ymesh,dose_zmesh,dose,handle);
            set(handle,'FaceColor',contourcolor(p),'EdgeColor','None','FaceAlpha',transl(p));
        end
        legend(legend_matrix);
        view(30,30);
        xlabel('X axis (cm)','Fontsize',12)
        ylabel('Y axis (cm)','Fontsize',12)
        zlabel('Z axis (cm)','Fontsize',12)
    end
    title(['Isosurface plot: ',inputname(1)],'interpreter','none','FontSize',16);
    % Plot VOIs
    dicomrt_rendervoi(VOI,voi2plot,1,1,1);
else
    error('dicomrt_isosurface: VOI or voi2plot is missing. Exit now !');
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function range = p12_range ( )

%*****************************************************************************80
%
%% P12_RANGE returns an interval bounding the root for problem 12.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real RANGE(2), the minimum and maximum values of
%    an interval containing the root.
%
  range(1) = - 4.0;
  range(2) =   4.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [out,varargout] = runningExtreme(input,nSamples,type)
% runningmin - Computes a running extreme of an input vector or matrix
% Optional file header info (to give more details about the function than in the H1 line)
% Syntax: out = runningmin(input,nSamples,type)
% input - vector or matrix of data to return the running min or max from
% nSamples - The size of the window to check for mins or maxes
% type - 'min','max', or 'both' (default) which type of extremes to return
% out - running minimum or maximum requested.  If both minimum is returned
% first, then maximum
%   Example
%  [runMin,runMax] = runningExtreme(data,31,'both')
%
% Subfunctions: vanherk
%   See also: min, max, sort
%% AUTHOR    : Dan Kominsky
%%

if nargin<2
	error('RunningExtreme:notEnoughInputs','Invalid Usage - must specify window size');
elseif nargin ==2
	type = 'both';
end

if ~mod(nSamples,2)
	warning('RunningExtreme:evenWindowSize','Running Min/Max is not meaningful for a windown with an even number of samples');
	nSamples = nSamples+1;
end

[input,dimShift]=shiftdim(input);

switch lower(type)
	case 'min'
		out=vanherk(input,nSamples,'min');% Create output for minimum
	case 'max'
		out=vanherk(input,nSamples,'max'); % Create output for maximum
	otherwise
		out=vanherk(input,nSamples,'max');% Create output for maximum
		varargout{1} = out; % Transfer to the second output
		out=vanherk(input,nSamples,'min');% Create output for minimum
end
if dimShift
	out = shiftdim(out,-dimShift);% if we transposed, undo to return the same shape as we got.
	if nargout>1
		varargout{1} = shiftdim(varargout{1},-dimShift);
	end
end

end

function Y = vanherk(X,N,TYPE)
%  VANHERK    Fast max/min 1D filter
%
%    Y = VANHERK(X,N,TYPE) performs the 1D max/min filtering of the row
%    vector X using a N-length filter.
%    The filtering type is defined by TYPE = 'max' or 'min'. This function
%    uses the van Herk algorithm for min/max filters that demands only 3
%    min/max calculations per element, independently of the filter size.
%
%    If X is a 2D matrix, each column will be filtered separately.
%    X can be uint8 or double. If X is uint8 the processing is quite faster, so
%    dont't use X as double, unless it is really necessary.
%

% Initialization
% if strcmp(direc,'col')
%    X = X';
% end
switch lower(TYPE)
	case 'max'
		maxfilt = 1;
	case 'min'
		maxfilt = 0;
	otherwise
		error([ 'TYPE must be ' char(39) 'max' char(39) ' or ' char(39) 'min' char(39) '.'])
end

% Correcting X size
fixsize = 0;
addel = 0;
if mod(size(X,1),N) ~= 0
	fixsize = 1;
	addel = N-mod(size(X,1),N);
	if maxfilt
		f = [X; -Inf*ones(addel,size(X,2)) ]; % Change from adding zeros
	else
		f = [X; Inf*ones([addel size(X,2)])];
		%       f = [X; repmat(X(end,:),addel,1)];  % Adds a replication of the end of the matrix
	end
else
	f = X;
end
lf = size(f,1); % # of elements in adjusted matrix
lx = size(X,1); % # of elements in original matrix
clear X

% Declaring aux. mat.
g = f;
h = g;

% Filling g & h (aux. mat.)
ig = (1:N:size(f,1)).'; % First element of each window
ih = ig + N - 1;    % Last element of each window

if maxfilt
	for i = 2 : N
		igold = ig;
		ihold = ih;

		ig = ig + 1;
		ih = ih - 1;

		g(ig,:) = max(f(ig,:),g(igold,:));
		h(ih,:) = max(f(ih,:),h(ihold,:));
	end
else
	for i = 2 : N
		igold = ig;
		ihold = ih;

		ig = ig + 1;
		ih = ih - 1;

		g(ig,:) = min(f(ig,:),g(igold,:));
		h(ih,:) = min(f(ih,:),h(ihold,:));
	end
end
clear f

if fixsize % If we had to pad the data
	if addel > (N-1)/2 % If the padding is more than half a zone
		ig =  (N : 1 : lf - addel + floor((N-1)/2)).' ;
		ih = ( 1 : 1 : lf-N+1 - addel + floor((N-1)/2)).';
		if maxfilt
			Y = [ g(1+ceil((N-1)/2):N-1,:);  max(g(ig,:), h(ih,:)) ];
		else
			Y = [ g(1+ceil((N-1)/2):N-1,:);  min(g(ig,:), h(ih,:)) ];
		end
	else
		ig = ( N : 1 : lf ).';
		ih = ( 1 : 1 : lf-N+1 ).';
		if maxfilt
			Y = [ g(1+ceil((N-1)/2):N-1,:);  max(g(ig,:), h(ih,:));  h(lf-N+2:lf-N+1+floor((N-1)/2)-addel,:) ];
		else
			Y = [ g(1+ceil((N-1)/2):N-1,:);  min(g(ig,:), h(ih,:));  h(lf-N+2:lf-N+1+floor((N-1)/2)-addel,:) ];
		end
	end
else % not fixsize (addel=0, lf=lx)
	ig = ( N : 1 : lx ).';
	ih = ( 1 : 1 : lx-N+1 ).';
	if maxfilt
		Y = [  g(N-ceil((N-1)/2):N-1,:); max( g(ig,:), h(ih,:) );  h(lx-N+2:lx-N+1+floor((N-1)/2),:) ];
	else
		Y = [  g(N-ceil((N-1)/2):N-1,:); min( g(ig,:), h(ih,:) );  h(lx-N+2:lx-N+1+floor((N-1)/2),:) ];
	end
end

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: Sable

repo: mcbench-benchmarks


text:
Eat = [E.Longitude.*111, E.Latitude.*111, E.Depth];
Eatlon = newt2.Longitude.*111;
Eatran = ran(:,1).*111;
Eattim = newt2.Date;
incr = a.Date;
eat = [((Eat(:,1).^2 + Eat(:,2).^2 + Eat(:,3)).^0.5)];
%
%
%
p = min(a.Date);

tim2 = [];
tim3 = [];

tim1 = [];
tic;
for p = min(a.Date):0.0000011:max(a.Date)

    match1 = find(a.Date< (p + 0.0000005) & a.Date> (p - 0.0000005));
    tim1 = [tim1; match1];

end
toc;

time = a(tim1,3);
Eatlon = [a(tim1,:)]; %[a(tim,1), a(tim,2), a(tim,7)];

for p = min(a.Date):0.001+0.01:max(a.Date)
    match2 = find(a.Date< (p + 0.0001) & a.Date> (p - 0.0001));
    tim2 = [tim2; match2];
end
Eat2 = [eat(tim2)];

for p = min(a.Date):0.001+0.02:max(a.Date)
    match3 = find(a.Date< (p + 0.0001) & a.Date> (p - 0.0001));
    tim3 = [tim3; match3];
end
Eat3 = [eat(tim3)];

Eat = [Eat1 Eat2 Eat3];

for m = 0:6

    time1(1:3202,(m+1)) = (decyear(E.Date)+ 0.005*m);
    %Eat(1:3202,(m)) = [Eat(:,1) eat(:,2)+0.005*m;
end

ea = find(decyear(E.Date)>min(decyear(E.Date))+0.05*4);
E = E(ea,:);

%
%
%

m = 1:5:2435; %size(pairdist,1);
figure;
%plot3(eat (m-2), eat(m-1,:), eat(m,:), 'k.', 'Markersize', 0.5);
%plot3(a((m+40),3), a(m+20,3), a(m,3), 'k.', 'Markersize', 0.5);
u = plot3(eat(m),eat(m+10),eat(m+20));%,'k.');
clear m;
%
%
%

k = 1:798;
u = [eat(k),eat(k+0.75),eat(k+1.5)];%,eat(k+2.25),eat(k+3),eat(k+3.75)];
figure;
%plot3(u(:,1),u(:,2), u(:,3), 'k.', 'Markersize', 0.5);%u(m-2,:)
xlabel('kt');
ylabel('kt+T');
zlabel('kt+2T');

k = 1:10:3012245;
figure;
plot3(pairdist(k),pairdist(k+20), pairdist(k+40));%,'k.', 'Markersize', 0.5);

%plot(E(n-1,1), E(n,1),'k.', 'Markersize', 0.5);
%plot(E(n-1,2), E(n,2),'k.', 'Markersize', 0.5);
%plot(E(n-1,7), E(n,7),'k.', 'Markersize', 0.5);

k = 1:1:3193;
E = [Eatlon(k),Eatlon(k+2),Eatlon(k+3), Eatlon(4), Eatlon(k+5), Eatlon(6), Eatlon(k+7), Eatlon(k+8)];

mEatlon = sum(Eatlon,1)/3202;
Eatlon = Eatlon-mEatlon;
HSig = figure;
plot(Eattim, Eatlon);
axis([2500 3000 -10 -15]);

[ps,freq] = spectrum(Eatlon,256,0,[],0.0634);
Hpws = figure;
plot(freq,ps(:,1));
axis([0 0.01 0 30]);
spectrum(Eatlon);

author: CelsoReyes


text:
%MIT IAP Radar Course 2011
%Resource: Build a Small Radar System Capable of Sensing Range, Doppler,
%and Synthetic Aperture Radar Imaging
%
%Gregory L. Charvat

%Process Range vs. Time Intensity (RTI) plot

%NOTE: set up-ramp sweep from 2-3.2V to stay within ISM band
%change fstart and fstop bellow when in ISM band

clear all;
close all;

%read the raw data .wave file here
[Y,FS,NBITS] = wavread('running_outside_20ms.wav');

%constants
c = 3E8; %(m/s) speed of light

%radar parameters
Tp = 20E-3; %(s) pulse time
N = Tp*FS; %# of samples per pulse
fstart = 2260E6; %(Hz) LFM start frequency for example
fstop = 2590E6; %(Hz) LFM stop frequency for example
%fstart = 2402E6; %(Hz) LFM start frequency for ISM band
%fstop = 2495E6; %(Hz) LFM stop frequency for ISM band
BW = fstop-fstart; %(Hz) transmti bandwidth
f = linspace(fstart, fstop, N/2); %instantaneous transmit frequency

%range resolution
rr = c/(2*BW);
max_range = rr*N/2;

%the input appears to be inverted
trig = -1*Y(:,1);
s = -1*Y(:,2);
clear Y;

%parse the data here by triggering off rising edge of sync pulse
count = 0;
thresh = 0;
start = (trig > thresh);
for ii = 100:(size(start,1)-N)
    if start(ii) == 1 & mean(start(ii-11:ii-1)) == 0
        %start2(ii) = 1;
        count = count + 1;
        sif(count,:) = s(ii:ii+N-1);
        time(count) = ii*1/FS;
    end
end
%check to see if triggering works
% plot(trig,'.b');
% hold on;si
% plot(start2,'.r');
% hold off;
% grid on;

%subtract the average
ave = mean(sif,1);
for ii = 1:size(sif,1);
    sif(ii,:) = sif(ii,:) - ave;
end

zpad = 8*N/2;

%RTI plot
figure(10);
v = dbv(ifft(sif,zpad,2));
S = v(:,1:size(v,2)/2);
m = max(max(v));
imagesc(linspace(0,max_range,zpad),time,S-m,[-80, 0]);
colorbar;
ylabel('time (s)');
xlabel('range (m)');
title('RTI without clutter rejection');

%2 pulse cancelor RTI plot
figure(20);
sif2 = sif(2:size(sif,1),:)-sif(1:size(sif,1)-1,:);
v = ifft(sif2,zpad,2);
S=v;
R = linspace(0,max_range,zpad);
for ii = 1:size(S,1)
    %S(ii,:) = S(ii,:).*R.^(3/2); %Optional: magnitude scale to range
end
S = dbv(S(:,1:size(v,2)/2));
m = max(max(S));
imagesc(R,time,S-m,[-80, 0]);
colorbar;
ylabel('time (s)');
xlabel('range (m)');
title('RTI with 2-pulse cancelor clutter rejection');

% %2 pulse mag only cancelor
% figure(30);
% clear v;
% for ii = 1:size(sif,1)-1
%     v1 = abs(ifft(sif(ii,:),zpad));
%     v2 = abs(ifft(sif(ii+1,:),zpad));
%     v(ii,:) = v2-v1;
% end
% S=v;
% R = linspace(0,max_range,zpad);
% for ii = 1:size(S,1)
%     S(ii,:) = S(ii,:).*R.^(3/2); %Optional: magnitude scale to range
% end
% S = dbv(S(:,1:size(v,2)/2));
% m = max(max(S));
% imagesc(R,time,S-m,[-20, 0]);
% colorbar;
% ylabel('time (s)');
% xlabel('range (m)');
% title('RTI with 2-pulse mag only cancelor clutter rejection');

author: lukeweston

repo: SimpleFMCWRadar


text:
function value = i4vec_any_negative ( n, a )

%*****************************************************************************80
%
%% I4VEC_ANY_NEGATIVE: ( any A < 0 ) for I4VEC's.
%
%  Discussion:
%
%    An I4VEC is a vector of I4's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 December 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries.
%
%    Input, integer A(N), the vector.
%
%    Output, logical VALUE is TRUE if any entry is negative.
%
  value = any ( a(1:n) < 0 );

  return
end


text:
% makegausslpfir        Gaussian LP filter
%
% win = makegausslpfir( Fc, Fs, s )
%
% Fc    corner frequency [Hz]
% Fs    sampling frequency [Hz]
% s     support [SD], minimum 3, default 4
%
% win   fir
%
% see also  firfilter

% 02-oct-13 ES

function gwin = makegausslpfir( Fc, Fs, s )

nargs = nargin;
if nargs < 1 || isempty( Fc ), Fc = 100; end
if nargs < 2 || isempty( Fs ), Fs = 1000; end
if nargs < 3 || isempty( s ), s = 4; end
s = max( s, 3 );

sd = Fs / ( 2 * pi * Fc );
x = -ceil( s * sd ) : ceil( s * sd );
gwin = 1/( 2 * pi * sd ) * exp( -( x.^2/2/sd.^2 ) );
gwin = gwin / sum( gwin );

return

% EOF

% to determine the proper support:
cwin = cumsum( gwin );
length( cwin ), [ find( cwin >= 0.001, 1, 'first' ) - 1, find( cwin <= 0.999, 1, 'last' ) + 1 ]
% 3 SD are usually enough

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function i4_factorial2_test ( )

%*****************************************************************************80
%
%% I4_FACTORIAL2_TEST tests I4_FACTORIAL2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 January 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4_FACTORIAL2_TEST\n' );
  fprintf ( 1, '  I4_FACTORIAL2 evaluates the double factorial function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '         N      Exact         I4_FACTORIAL2(N)\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, n, f1 ] = i4_factorial2_values ( n_data );

    if ( n_data == 0 )
      break
    end

    f2 = i4_factorial2 ( n );

    fprintf ( 1, '  %8d  %12d  %12d\n', n, f1, f2 );

  end

  return
end


text:
function m = marginal_family(engine, n, add_ev)
% MARGINAL_FAMILY Compute the marginal on i's family (loopy)
% m = marginal_family(engine, n, add_ev)

if nargin < 3, add_ev = 0; end

bnet = bnet_from_engine(engine);
ns = bnet.node_sizes;
ps = parents(bnet.dag, n);
dom = [ps n];
CPD = bnet.CPD{bnet.equiv_class(n)};

switch engine.msg_type
  case 'd',
   % The method is similar to the following HMM equation:
   % xi(i,j,t) = normalise( alpha(i,t) * transmat(i,j) * obsmat(j,t+1) * beta(j,t+1) )
   % where xi(i,j,t) = Pr(Q(t)=i, Q(t+1)=j | y(1:T))
   % beta == lambda, alpha == pi, alpha from each parent = pi msg
   % In general, if A,B are parents of C,
   % P(A,B,C) = P(C|A,B) pi_msg(A->C) pi_msg(B->C) lambda(C)
   % where lambda(C) = P(ev below and including C|C) = prod incoming lamba_msg(children->C)
   % and pi_msg(X->C) = P(X|ev above) etc

   T = dpot(dom, ns(dom), CPD_to_CPT(CPD));
   for j=1:length(ps)
     p = ps(j);
     pi_msg = dpot(p, ns(p), engine.msg{n}.pi_from_parent{j});
     T = multiply_by_pot(T, pi_msg);
   end
   lambda = dpot(n, ns(n), engine.msg{n}.lambda);
   T = multiply_by_pot(T, lambda);
   T = normalize_pot(T);
   m = pot_to_marginal(T);
   if ~add_ev
     m.T = shrink_obs_dims_in_table(m.T, dom, engine.evidence);
   end
 case 'g',
  if engine.disconnected_nodes_bitv(n)
    m.T = 1;
    m.domain = dom;
    if add_ev
      m = add_ev_to_dmarginal(m, engine.evidence, ns)
    end
    return;
  end

  [m, C, W] = gaussian_CPD_params_given_dps(CPD, dom, engine.evidence);
  cdom = myintersect(dom, bnet.cnodes);
  pot = linear_gaussian_to_cpot(m, C, W, dom, ns, cdom, engine.evidence);
  % linear_gaussian_to_cpot will set the effective size of observed nodes to 0,
  % so we need to do this explicitely for the messages, too,
  % so they are all the same size.
  obs_bitv = ~isemptycell(engine.evidence);
  ps = parents(engine.msg_dag, n);
  for j=1:length(ps)
    p = ps(j);
    msg = engine.msg{n}.pi_from_parent{j};
    if obs_bitv(p)
      pi_msg = mpot(p, 0);
    else
      pi_msg = mpot(p, ns(p), 0, msg.mu, msg.Sigma);
    end
    pot = multiply_by_pot(pot, mpot_to_cpot(pi_msg));
  end
  msg = engine.msg{n}.lambda;
  if obs_bitv(n)
    lambda = cpot(n, 0);
  else
    lambda = cpot(n, ns(n), 0, msg.info_state, msg.precision);
  end
  pot = multiply_by_pot(pot, lambda);
  m = pot_to_marginal(pot);
  if add_ev
    m = add_evidence_to_gmarginal(m, engine.evidence, bnet.node_sizes, bnet.cnodes);
  end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [status] = putexif(dat,fname, refname)
%[status] = putexif(dat,fname, refname) save and image with Exif data from
%                  a reference image file
% dat     = array of image data
% fname   = name of image file to be saved, e.g., 'myfile.jpg'
% refname = name (path) of reference image file
%
%Needs ExifTool.exe, written by Phil Harvey
% http://www.sno.phy.queensu.ca/~phil/exiftool/
% 1. Download the file, exiftool(-k).exe
% 2. Rename this file, exiftool.exe
% 3. Save this file in a folder on your Matlab path, e.g. .../matlab
%Peter Burns, 28 May 2013
%             22 July 2013, following suggestions from jhh and Jonathan via
%                           Matlab Central.

test = which('exiftool.exe');
if isempty(test)
    disp('ExifTool not available:');
    disp('Please download from,')
    disp('   http://www.sno.phy.queensu.ca/~phil/exiftool/')
    disp('or make sure that the installed exiftool.exe is on your Matlab path')
    beep
    status =1;
    return
else

    [exifdata, nf] = getexif(refname);
   % Save data as image file with Matlab metadata
    imwrite(dat, fname, 'quality', 100);
    [exifdata2, nf2] = getexif(fname);
   % Replace matadata from this file with desired tags from reference file
   % temp1=['exiftool -m -tagsfromfile ',refname,' -all:all ', fname];
    temp1=['"' test '" -m -tagsfromfile "',refname,'" -all:all "', fname, '"'];

    %[status, junk] = system(temp1);
    %Second output variable appears to suppress '1 image files updated'
    % being returned to Command window (~ also works).
    [status, junk] = system(temp1);

   % Approximate check
    [exifdata3, nf3] = getexif(fname);

   % Delete extra copy of reference image file
    temp3 = ['del "',fname,'_original"'];
    status = system(temp3);

   % Test approx. number of tags
    if abs(nf3-nf2<10)
        disp('Warning: Exif tags may not have been copied');
        status = 1;
    end

end


text:
function cx = csscal ( n, sa, cx, incx )

%*****************************************************************************80
%
%% CSSCAL scales a complex vector by a real constant.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979,
%    ISBN13: 978-0-898711-72-1,
%    LC: QA214.L56.
%
%    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
%    Basic Linear Algebra Subprograms for Fortran Usage,
%    Algorithm 539,
%    ACM Transactions on Mathematical Software,
%    Volume 5, Number 3, September 1979, pages 308-323.
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, real SA, the multiplier.
%
%    Input, complex CX(*), the vector to be scaled.
%
%    Input, integer INCX, the increment between successive entries of
%    the vector CX.
%
%    Output, complex CX(*), the scaled vector.
%
  cx(1:incx:1+(n-1)*incx) = sa * cx(1:incx:1+(n-1)*incx);

  return
end


text:
function y = hmean(x)
% HMEAN   1/mean(1/input)
%
%   Description:
%   Y = HMEAN(X) evaluates y=1./mean(1./x)

% Copyright (c) 1998 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

y=1./mean(1./x);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function  [bdLat,bdLng]= gcj2bd(gcjLat, gcjLng)

    bdLat = gcjLat;
    bdLng = gcjLng;

    x_pi = pi * 3000.0 / 180.0;
    inChina = ~outOfChina(gcjLat, gcjLng);
    if ~any(inChina),return;end

    x = gcjLng(inChina);
    y = gcjLat(inChina);
    z = hypot(x, y) + 0.00002 * sin(y * x_pi);
    theta = atan2(y, x) + 0.000003 * cos(x * x_pi);
    bdLng(inChina) = z.* cos(theta) + 0.0065;
    bdLat(inChina) = z.* sin(theta) + 0.006;

end

author: googollee

repo: eviltransform


text:
function K = white_noise_correlation ( s, t )

%*****************************************************************************80
%
%% WHITE_NOISE_CORRELATION evaluates the white_noise correlation function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Petter Abrahamsen,
%    A Review of Gaussian Random Fields and Correlation Functions,
%    Norwegian Computing Center, 1997.
%
%  Parameters:
%
%    Input, real S(*), T(*), pairs of argument values.
%
%    Output, real K(*), the correlation function values
%
  K = zeros ( size ( s ) );

  i = find ( s == t );

  K(i) = 1.0;

  return
end


text:
function C = xeucn( A, T, shape )
% n-dimensional euclidean distance between each window in A and template T.
%
% Similar to normxcorrn, except at each point (i,j) calculates the
% euclidean distance between the T and the window in A surrounding the
% point, storing the result in C(i,j).
%
% USAGE
%  C = xeucn( A, T, [shape] )
%
% INPUTS
%  A           - first d-dimensional matrix
%  T           - second d-dimensional matrix
%  shape       - ['full'] 'valid', or 'same' (see convn)
%
% OUTPUTS
%  C           - correlation matrix
%
% EXAMPLE
%  T=gaussSmooth(rand(20),2); A=repmat(T,[3 3]);
%  C1=normxcorrn(T,A);  C2=xcorrn(A,T);  C3=xeucn(A,T);
%  figure(1); im(C1);  figure(2); im(C2);  figure(3); im(-C3);
%
% See also XCORRN, CONVNFAST
%
% Piotr's Image&Video Toolbox      Version 2.0
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

if( nargin<3 || isempty(shape)); shape='full'; end
nd = ndims(T);
if(nd~=ndims(A)); error('xeucn: T and A must have same ndims'); end;

% flip for conv purposes
if(nd==2); T=rot90(T,2); else for d=1:nd; T=flipdim(T,d); end; end

% The expression for euclidean distance can be rewritten as:
%   D(k,l) = sumj( (Akj - Tlj).^2 )
%          = sumj( Akj.^2 ) + sumj( Tlj.^2 ) - 2*sumj(Akj.*Tlj);
% T is constant.  Hence simply need square of A in each window, as
% well as each dot product between A and T.
Amag = localSum( A.*A, size(T), shape );     % sum of squares per window
Tmag = T.^2;  Tmag = sum( Tmag(:) );         % sum of squares of T
C = Amag + Tmag - 2 * convnFast(A,T,shape);  % distance squared
% C( Amag<.01 ) = Tmag;  % prevent numerical errors
C = real(sqrt(real(C)));

author: zouchuhang

repo: LayoutNet


text:
function cb = buildCFun(mode,sobj,svar,opts,ncon)
% BUILDCFUN Build a C Code function for nonlinear callbacks

if(nargin > 3 && ~isempty(opts) && isfield(opts,'cbmode') && strcmpi(opts.cbmode,'cppad'))
    tname = 'symb_cadtemp.cpp';
    fname = 'symb_ccb.cpp';
    cmode = 'A';
else
    tname = 'symb_ctemp.c';
    fname = 'symb_ccb.c';
    cmode = 'C';
end
cb = [];

%Determine callback type
switch(mode)
    case 'new'
        %If compiling against CppAD, must use Intel or Visual Studio compiler
        if(cmode=='A')
            SymBuilder.CheckCppADCompile(true);
        elseif(cmode=='C') %check for LCC, doesn't work (don't know why)
            c = mex.getCompilerConfigurations('c');
            if(~isempty(c) && ~isempty(strfind(c.Name,'Lcc')))
                error('Lcc does not correctly compile SymBuilder C-Code callbacks. Please use Visual Studio or Windows SDK.');
            end
        end
        %Find template
        p = which(tname);
        if(isempty(p)), error('Cannot find C template'); end
        %Copy template
        str = fileread(p);
        s2 = sprintf('Symbolic Builder Auto-Generated Callback Function\n');
        s2 = sprintf('%s * Generated %s',s2,datestr(now));
        str = regexprep(str,'SYMB_CTEMP - Template for generating SymBuilder C Code Callbacks',s2);
        fp = fopen(fname,'w');
        fprintf(fp,'%s\n\n',str);
        fclose(fp);
        return;
    case 'obj'
        title = 'Objective Function Callback';
        if(cmode=='A')
            fcall = sprintf('template <typename Type> Type objective(const vector<Type> &x)');
        else
            fcall = sprintf('double objective(double *x)');
        end
        var = 'j';
        cb = @(x) symb_ccb('obj',x);
        if(opts.verbose), fprintf('Generating Objective....'); end
    case 'grad'
        title = 'Objective Gradient Callback';
        fcall = sprintf('void gradient(double *x, double *v)');
        var = 'g';
        cb = @(x) symb_ccb('grad',x);
        if(nargin > 1 && ~isempty(sobj))
            if(opts.verbose), fprintf('Generating Gradient....'); end
        end
    case 'con'
        title = 'Constraint Function Callback';
        if(cmode=='A')
            fcall = sprintf('template <typename Type> void constraints(const vector<Type> &x, vector<Type> &v)');
        else
            fcall = sprintf('void constraints(double *x, double *v)');
        end
        var = 'c';
        cb = @(x) symb_ccb('con',x);
        if(nargin > 1 && ~isempty(sobj))
            if(opts.verbose), fprintf('Generating Constraints....'); end
        end
    case 'jac'
        title = 'Constraint Jacobian Callback';
        fcall = sprintf('void jacobian(double *x, double *pr, mwIndex *ir, mwIndex *jc)');
        var = 'J';
        cb = @(x) symb_ccb('jac',x);
        if(nargin > 1 && ~isempty(sobj))
            if(opts.verbose), fprintf('Generating Jacobian....'); end
        end
    case 'hess'
        title = 'Hessian of the Lagrangian Callback';
        fcall = sprintf('void hessian(double *x, double sigma, double *lambda, double *pr, mwIndex *ir, mwIndex *jc)');
        var = 'H';
        cb = @(x,sigma,lambda) symb_ccb('hess',x,sigma,lambda);
       if(nargin > 1 && ~isempty(sobj))
           if(opts.verbose), fprintf('Generating Hessian....'); end
       end
    otherwise
        error('Unknown MFun mode: %s',mode);
end

%Open file
fp = fopen(fname,'a+');
%Write Header
fprintf(fp,'// %s\n',title);
fprintf(fp,fcall);
fprintf(fp,'\n{\n');

%If we have a function to write
if(nargin > 1 && ~isempty(sobj))
    %Need to rename all variables to x[0]...x[n]
    xvar = cell(length(svar),1);
    for i = 1:length(xvar)
        xvar{i} = sprintf('x[%d]',i);
    end
    %Ensure both columns
    if(size(svar,1) > 1), svar = svar.'; end
    if(size(xvar,1) > 1), xvar = xvar'; end
    %Subs out individual symbolic variables into our indexed list and converts to normal numbers
    wstate = warning('off','symbolic:sym:sym:DeprecateExpressions');
    eq = vpa(subs(sobj,svar,xvar),16); %this takes too long - any suggestions?
    warning(wstate);
    %Sub out Lambda if Hessian
    if(var=='H')
        l = cell(ncon,1);
        l2 = cell(ncon,1);
        for i = 1:ncon
            l{i} = sprintf('lambda(%d)',i);
            l2{i} = sprintf('lambda[%d]',i);
        end
        wstate = warning('off','symbolic:sym:sym:DeprecateExpressions');
        eq = subs(eq,l,l2);
        warning(wstate);
    end

    %Enter Equations (Var Type Dictates Entry Type)
    switch(var)
        case 'j' %scalar
            fprintf(fp,'%s\n',ppc(cmode,regexprep(ccode(eq),'t0 = ','return ')));
            %Spaces at end
            fprintf(fp,'}\n\n');
            %Extra #Var function
            writeVarFcn(fp,length(svar));
        case {'g','c'} %vector, dense
            %HACK
            if(~isa(eq,'sym')), eq = sym(eq); end
            if(size(eq,2) > 1), eq = eq.'; end
%             %Convert equations to string, write each as we go
%             for i = 1:length(eq)
%                 fprintf(fp,'%s\n',ppc(cmode,regexprep(ccode(eq(i)),'t0 = ',sprintf('v[%d] = ',i-1))));
%             end
            if(length(eq) > 1) %write all equations at once
                fprintf(fp,'%s\n',ppc(cmode,regexprep(ccode(eq),'eq[(\d*)\][0\]','v[$1\]')));
            else
                fprintf(fp,'%s\n',ppc(cmode,regexprep(ccode(eq),'t0 = ','v[0] = ')));
            end
            if(var=='c')
                %Spaces at end
                fprintf(fp,'}\n\n');
                %Extra #Con function
                writeConFcn(fp,length(eq));
            end
        case {'J','H'} %matrix, sparse
            %Get NonZero elements from equation
            nzel = logical(eq ~= 0);
            [rows,cols] = find(sparse(nzel));
            nz = nnz(nzel);
            nzsym = eq(nzel);
            %Write Ir (rows)
            for i = 1:length(rows)
                fprintf(fp,'  ir[%d] = %d;\n',i-1,rows(i)-1);
            end
            %Write Jc (col starts)
            s = 0;
            for i = 1:size(eq,2)+1
                fprintf(fp,'  jc[%d] = %d;\n',i-1,s);
                s = s + sum(cols==i);
            end
            %Write Pr (vals/eqs)
%             for i = 1:length(nzsym)
%                 fprintf(fp,'%s\n',regexprep(ccode(nzsym(i)),'t0 = ',sprintf('pr[%d] = ',i-1)));
%             end
            if(length(nzsym) > 1)
                str = regexprep(ccode(nzsym),'nzsym[(\d*)\][0\]','pr[$1\]');
                if ~isempty(strfind(str, 'nzsym'))
                    nzsym = nzsym.';
                    str = regexprep(ccode(nzsym),'nzsym[(\d*)\][0\]','pr[$1\]');
                end
                fprintf(fp,'%s\n',ppc(cmode,str));
            else
                fprintf(fp,'%s\n',regexprep(ccode(nzsym),'t0 = ','pr[0] = '));
            end
            %Write extra NNZ function
            fprintf(fp,'}\n\n');
            writeNNZFcn(fp,nz,mode);
    end
else
    fprintf(fp,'   mexErrMsgTxt("%s is not supported in this SymBuilder Callback");\n',mode);
end
%Spaces at end
fprintf(fp,'}\n\n');
%NNZ default functions
if((nargin < 2 || isempty(sobj)) && (strcmpi(mode,'jac') || strcmpi(mode,'hess')))
    writeNNZFcn(fp,0,mode);
    fprintf(fp,'}\n\n');
end
%Default con function
if((nargin < 2 || isempty(sobj)) && strcmpi(mode,'con'))
    writeConFcn(fp,0);
    fprintf(fp,'}\n\n');
end
%Close file
fclose(fp);
if(nargin > 1 && ~isempty(sobj))
    if(opts.verbose), fprintf('Done\n'); end
end

%Write getNoVarFcn
function writeVarFcn(fp,len)
fprintf(fp,'//Return number of variables\n');
fprintf(fp,'mwIndex getNoVar()\n{\n');
fprintf(fp,'   return %d;\n',len);

%Write getNoConFcn
function writeConFcn(fp,len)
fprintf(fp,'//Return number of constraints\n');
fprintf(fp,'mwIndex getNoCon()\n{\n');
fprintf(fp,'   return %d;\n',len);

%Write getNNZFcn
function writeNNZFcn(fp,nz,mode)
if(strcmpi(mode,'jac'))
    fprintf(fp,'//Return number of nonzeros in the Jacobian\n');
    fprintf(fp,'mwIndex getNNZJac()\n{\n');
else
    fprintf(fp,'//Return number of nonzeros in the Hessian\n');
    fprintf(fp,'mwIndex getNNZHess()\n{\n');
end
fprintf(fp,'   return %d;\n',nz);

%Preprocess String to make it compatible with CppAD
function str = ppc(mode,str)

if(strcmpi(mode,'C')), return; end
if(isempty(strfind(str,'.0'))), return; end

str = strrep(str,'.0)',')');
str = strrep(str,'.0,',',');
str = strrep(str,'.0;',';');
str = strrep(str,'.0+','+');
str = strrep(str,'.0-','-');
str = strrep(str,'.0*','*');
str = strrep(str,'.0/','/');
str = strrep(str,'.0 ',' ');

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [c, ceq] = ma_radiusConstraint(stateLog, eventID, lbR, ubR, bodyIDApply, celBodyData, maData)
%ma_radiusConstraint Summary of this function goes here
%   Detailed explanation goes here
    normFact = 10;

    if(ischar(eventID) && strcmpi(eventID,'final'))
        eventNum = max(stateLog(:,13));
    else
%         hMAMainGUI = findall(0,'tag','ma_MainGUI');
%         maData = getappdata(hMAMainGUI,'ma_data');
        [~, eventNum] = getEventByID(eventID, maData.script);
    end

    eventLog = stateLog(stateLog(:,13)==eventNum,:);

    bodyEventLog = eventLog(eventLog(:,8)==bodyIDApply,:);
    if(isempty(bodyEventLog))
        finalEntry = eventLog(end,:);
    else
        finalEntry = bodyEventLog(end,:);
    end

    bodyID = finalEntry(8);

    if(bodyID == bodyIDApply || bodyIDApply==-1)
        rVect = finalEntry(2:4)';
        r = norm(rVect);

        if(lbR == ubR)
            c = [0 0];
            ceq(1) = r - ubR;
        else
            c(1) = lbR - r;
            c(2) = r - ubR;
            ceq = [0];
        end
        c = c/normFact;
        ceq = ceq/normFact;
    else
        c = [0 0];
        ceq = [0];
    end

end


text:
clc
clear all
tic
%***********************User inserted input for the problem***************

nosp = 5;             %reads the number of sales points considered
noft = 4;                %reads the number of factories
nopd = 3;         %number of planning period
trkcap = 70;                    %Capacity of the truck

%***********************inputs from excel file for the formulation ***************

demand = [35 35	60	80	100; 25	100	125	55	20; 150	120	100	60	110; 50	70	40	140	40 ; 90	110	90	50	30;0	110	90	40	50; 60	70	0	110	10; 140	70	90	100	30 ;10	130	130	40	140; 30	90	70	40	0; 70	100	100	50	90; 110	50	50	110	100 ]; %stores the demand in the time horizons for each sales point
lostdemandcost = [30	35	25	40	55] ;  %lost demand per unit for each sales point
transcost = [ 100	400	800	1100	1700 ; 1500	1800	1200	2200	1500 ; 1400	1000	400	1000	700 ; 800	700	900	700	1200] ; %xlsread('greedyinput','Sheet1','O3:S6');        %transportation cost matrix
setupcost = [ 4500	2000	2500	3000 ] ;  %setup costs for individual factories
factorycap = [300	500	250	150]; % xlsread('greedyinput','Sheet2','B4:E4');       %capacity of the factories
factholdcost = [5	2	4	3]; %xlsread('greedyinput','Sheet2','B5:E5');     %holding cost per unit in a factory
for i = 1:nopd
quantityshipped{i} = zeros(noft,nosp);
factoryresiduals{i} = zeros(1,noft);
factoryactivation{i} = zeros(1,noft);
salespointresiduals{i} = zeros(1,nosp);
end

%*****************************counter-check if the inputs are well read*********

disp(demand);
disp(lostdemandcost);
disp(transcost);
disp(setupcost);
disp(factorycap);
disp(factholdcost);
for i = 1:nopd
disp(quantityshipped{i});
end

%****************************starts the calculation of the greedy alg********

initlostsales = calclostsales(demand,lostdemandcost,nopd,nosp); %Calculates the lost demand with the initial assumption that all the demand is lost
disp(initlostsales);
cumsetupcost = 0;
cumtransportcost = 0;
cumholdingcost = 0;
saveddemand = 0;
salespoint{nosp}= [0,0];
for ft = 1:noft                                          %Initialize the factory activation status and residual for each time preiod (zero for the first one)
           factory{ft}(1,1) = 0;
           factory{ft}(1,2) = 0;
end
for time = 1:nopd                                               %iterates for each planning period beginning from the 1st

       for sp = 1:nosp
           salespoint{sp}(1,1) = 0;                             %the array stores the activation status and the demand of a sales point
            salespoint{sp}(1,2) =  demand(time,sp);             %this populates the sales point for a time horizon with demand and filling status
       end

       for ft = 1:noft                                          %Initialize the factory activation status and residual for each time preiod (zero for the first one)
           factory{ft}(1,1) = 0;

       end

   % loadfactory(capacity,prevresidual);

    for points = 1:nosp
        [selectedsp,maxlostsales] = choosesalespt(lostdemandcost,salespoint,nosp);             %Chooses the salespoint with the highest possible lost sales
        disp('----------------------------');
        fprintf('Selected Sales Point: %d\n', selectedsp);
        fprintf('Potential Lost Sales cost: %.2f\n', maxlostsales);

        salespoint{selectedsp}(1,1)=1;

            [costunit,selectedf] = choosefactory(factorycap,factholdcost,transcost,setupcost,factory,trkcap,salespoint,selectedsp); %chooses the preferable factory to deliver to the selected salespoint
            fprintf('selected factory: %d\n',selectedf);
            fprintf('Cost of the selected factory: %.2f\n', costunit);
            quantityshipped{time}(selectedf,selectedsp) = salespoint{selectedsp}(1,2);
            if (factory{selectedf}(1,2)>= salespoint{selectedsp}(1,2))
            factory{selectedf}(1,2) = factory{selectedf}(1,2) - salespoint{selectedsp}(1,2);

            % saveddemand = saveddemand + salespoint{selectedsp}(1,2)*lostdemandcost(selectedsp);                             % accumulates the reduction of the lost demand for each sent item to a destination sales point
            end
            if(factory{selectedf}(1,1)==0 && factory{selectedf}(1,2)<= salespoint{selectedsp}(1,2))
                factory{selectedf}(1,2) = factory{selectedf}(1,2)+factorycap(selectedf)-salespoint{selectedsp}(1,2);
                cumsetupcost = cumsetupcost + setupcost(selectedf);                                                         %accumulates the activation cost for each factory when there is an accompanying activation
               % cumtransportcost = cumtransportcost + ceil(salespoint{selectedsp}(1,2)/trkcap)*transcost(selectedf,selectedsp);
                factory{selectedf}(1,1)=1;
            end
            saveddemand = saveddemand + salespoint{selectedsp}(1,2)*lostdemandcost(selectedsp);
            cumtransportcost = cumtransportcost + ceil(salespoint{selectedsp}(1,2)/trkcap)*transcost(selectedf,selectedsp); %accumulates the transportaton cost for each iteration

            fprintf('Cumulative setupcost: %.2f\n',cumsetupcost);
            fprintf('Cumulative trasportation cost: %.2f\n',cumtransportcost);
    end
    currentholdingcost = 0;
    for i = 1:noft
        currentholdingcost= currentholdingcost + factory{i}(1,2)*factholdcost(i);                                         %the holding costs are accumulatedat each factory at the end of each planning period
        factoryresiduals{time}(i) = factory{i}(1,2);
        factoryactivation{time}(i) = factory{i}(1,1);
        factory{i}(1,1)= 0;
    end
    cumholdingcost = cumholdingcost + currentholdingcost;                                                                  % holding cost of the whole factory at the end of a planning period
    objectivesolution = initlostsales -saveddemand+ cumholdingcost+cumsetupcost+cumtransportcost
    end
%objectivesolution = initlostsales-saveddemand+cumsetupcost+cumtransportcost
%objectivesolution = initlostsales -saveddemand+ cumholdingcost+cumsetupcost+cumtransportcost
disp('====================shipped quantity matrix==================')
for i = 1:nopd
disp(quantityshipped{i});
end
disp('=====================factory residual matrix=================')
for i = 1:nopd
disp(factoryresiduals{i});
end
disp('=====================factory activation status===============')
for i = 1:nopd
disp(factoryactivation{i});
end
objectivesolution
% New_transportation_costs=0;
% F_holding_cost=0;
% Act_cost= 0;
% for t=1:nopd
%     for i=1:noft
%         F_holding_cost = F_holding_cost + factholdcost(i) *factoryresiduals{t}(1,i);
%         Act_cost = Act_cost + factoryactivation{t}(1,i)* setupcost(i);
%
%
%         for j =1:nosp
%     New_transportation_costs= ceil(quantityshipped{t}(i,j)/trkcap) * transcost(i,j)+ New_transportation_costs;
%         end
%     end
% end
% New_transportation_costs
% F_holding_cost
% Act_cost
% tot_cost = New_transportation_costs + F_holding_cost +Act_cost
% saveddemand
% initlostsales
% cumtransportcost
toc


text:
%SerialLink.IKINEM Numerical inverse kinematics by minimization
%
% Q = R.ikinem(T) is the joint coordinates corresponding to the robot
% end-effector pose T which is a homogenenous transform.
%
% Q = R.ikinem(T, Q0, OPTIONS) specifies the initial estimate of the joint
% coordinates.
%
% In all cases if T is 4x4xM it is taken as a homogeneous transform sequence
% and R.ikinem() returns the joint coordinates corresponding to each of the
% transforms in the sequence.  Q is MxN where N is the number of robot joints.
% The initial estimate of Q for each time step is taken as the solution
% from the previous time step.
%
% Options::
% 'pweight',P      weighting on position error norm compared to rotation
%                  error (default 1)
% 'stiffness',S    Stiffness used to impose a smoothness contraint on joint
%                  angles, useful when N is large (default 0)
% 'qlimits'        Enforce joint limits
% 'ilimit',L       Iteration limit (default 1000)
% 'nolm'           Disable Levenberg-Marquadt
%
% Notes::
% - PROTOTYPE CODE UNDER DEVELOPMENT, intended to do numerical inverse kinematics
%   with joint limits
% - The inverse kinematic solution is generally not unique, and
%   depends on the initial guess Q0 (defaults to 0).
% - The function to be minimized is highly nonlinear and the solution is
%   often trapped in a local minimum, adjust Q0 if this happens.
% - The default value of Q0 is zero which is a poor choice for most
%   manipulators (eg. puma560, twolink) since it corresponds to a kinematic
%   singularity.
% - Such a solution is completely general, though much less efficient
%   than specific inverse kinematic solutions derived symbolically, like
%   ikine6s or ikine3.% - Uses Levenberg-Marquadt minimizer LMFsolve if it can be found,
%   if 'nolm' is not given, and 'qlimits' false
% - The error function to be minimized is computed on the norm of the error
%   between current and desired tool pose.  This norm is computed from distances
%   and angles and 'pweight' can be used to scale the position error norm to
%   be congruent with rotation error norm.
% - This approach allows a solution to obtained at a singularity, but
%   the joint angles within the null space are arbitrarily assigned.
% - Joint offsets, if defined, are added to the inverse kinematics to
%   generate Q.
% - Joint limits become explicit contraints if 'qlimits' is set.
%
% See also fminsearch, fmincon, SerialLink.fkine, SerialLink.ikine, tr2angvec.

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

function qt = ikinem(robot, tr, varargin)

    opt.pweight = 1;
    opt.stiffness = 0;
    opt.qlimits = false;
    opt.ilimit = 1000;
    opt.lm = true;
    opt.col = 2;

    [opt,args] = tb_optparse(opt, varargin);

    % check if optional argument is a valid q
    q0 = args{1};
    if numel(q0) ~= robot.n
        error('q0 length must match number of joints in robot');
    end

    for i=1:size(tr,3)
        T = tr(:,:,i);

        if opt.qlimits
            % constrained optimization to handle joint limits
            options = optimset('MaxIter', opt.ilimit);
            qlim = robot.qlim;

            [q, ef, exflag, output] = fmincon( @(x) costfun(x, robot, T, opt), q0, ...
                [], [], [], [], ...
                qlim(:,1), qlim(:,2), ...
                [], options);

            if opt.verbose
                fprintf('final error %f, %d iterations, %d evalations\n', ...
                    ef, output.iterations, output.funcCount);
            end

        else
            % no joint limits, unconstrained optimization
            if exist('LMFsolve') == 2 && opt.lm
                [q, ef, count] = LMFsolve( @(x) costfun(x, robot, T, opt), q0, 'MaxIter', opt.ilimit);
                q = q';
                if opt.verbose
                    fprintf('final error %f, %d iterations\n', ...
                        ef, count);
                end
            else
                options = optimset('MaxIter', opt.ilimit);
                [q, ef, exflag, output] = fminsearch( @(x) costfun(x, robot, T, opt), q0, options);

                if opt.verbose
                    fprintf('final error %f, %d iterations, %d evalations\n', ...
                        ef, output.iterations, output.funcCount);
                end
            end
        end

        qt(i,:) = q;
    end

    if opt.verbose
        robot.fkine(qt)
    end
end

% The cost function, this is the value to be minimized
function E = costfun(q, robot, T, opt)

    Tq = robot.fkine(q);
    % find the pose error in SE(3)
    dT = transl(T) - transl(Tq);

    % translation error
    E = norm(dT) * opt.pweight;

    % rotation error
    %  find dot product of
    dd = dot(T(1:3,opt.col), Tq(1:3,opt.col));
    %E = E + (1 - dd)^2*100000 ;
    E = E + acos(dd)^2*1000 ;

    if opt.stiffness > 0
        % enforce a continuity constraint on joints, minimum bend
        E = E + sum( diff(q).^2 ) * opt.stiffness;
    end
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function varargout = ReadData3D(varargin)
% This function ReadData3D allows the user to open medical 3D files. It
% supports the following formats :
%
%   Dicom Files ( .dcm , .dicom )
%   V3D Philips Scanner ( .v3d )
%   GIPL Guys Image Processing Lab ( .gipl )
%   HDR/IMG Analyze ( .hdr )
%   ISI Files ( .isi )
%   NifTi ( .nii )
%   RAW files ( .raw , .* )
%   VMP BrainVoyager ( .vmp )
%   XIF HDllab/ATL ultrasound ( .xif )
%   VTK Visualization Toolkit ( .vtk )
%   Insight Meta-Image ( .mha, .mhd )
%   Micro CT ( .vff )
%   PAR/REC Philips ( .par, .rec)
%
% usage:
%
% [V,info]=ReadData3D;
%
% or,
%
% [V,info]=ReadData3D(filename)
%
% or,
%
% [V,info]=ReadData3D(filename,real);
%
%
% outputs,
%   V : The 3D Volume
%   info : Struct with info about the data
%        Always the following fields are present
%        info.Filename : Name of file
%        info.Dimensions : Dimensions of Volume
%        info.PixelDimensions : Size of one pixel / voxel
%   real : If set to true (default), convert the raw data to
%        type Single-precision  and rescale data to real units
%        (in CT Hounsfield). When false, it returns the raw-data.
%
% Warning!
%  The read functions are not fully implemented as defined in
%  the file-format standards. thus do not use this function for
%  critical applications.
%
%
% Function is written by D.Kroon University of Twente (July 2010)

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

% Last Modified by GUIDE v2.5 09-Nov-2010 14:12:50

% Begin initialization code - DO NOT EDIT
gui_Singleton = 0;
gui_State = struct('gui_Name',       mfilename, ...
    'gui_Singleton',  gui_Singleton, ...
    'gui_OpeningFcn', @ReadData3D_OpeningFcn, ...
    'gui_OutputFcn',  @ReadData3D_OutputFcn, ...
    'gui_LayoutFcn',  [] , ...
    'gui_Callback',   []);
if (nargin>2) && 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 ReadData3D is made visible.
function ReadData3D_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 ReadData3D (see VARARGIN)

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

% Update handles structure
guidata(hObject, handles);

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

%---- Start supported file formats ----%
data.fileformat(1).ext='*.dcm';
data.fileformat(1).type='Dicom Files';
data.fileformat(1).folder='dicom';
data.fileformat(1).functioninfo='dicom_read_header';
data.fileformat(1).functionread='dicom_read_volume';

data.fileformat(2).ext='*.gipl';
data.fileformat(2).type='GIPL Guys Image Processing Lab';
data.fileformat(2).folder='gipl';
data.fileformat(2).functioninfo='gipl_read_header';
data.fileformat(2).functionread='gipl_read_volume';

data.fileformat(3).ext='*.hdr';
data.fileformat(3).type='HDR/IMG Analyze';
data.fileformat(3).folder='hdr';
data.fileformat(3).functioninfo='hdr_read_header';
data.fileformat(3).functionread='hdr_read_volume';

data.fileformat(4).ext='*.isi';
data.fileformat(4).type='ISI Files';
data.fileformat(4).folder='isi';
data.fileformat(4).functioninfo='isi_read_header';
data.fileformat(4).functionread='isi_read_volume';

data.fileformat(5).ext='*.nii';
data.fileformat(5).type='NifTi';
data.fileformat(5).folder='nii';
data.fileformat(5).functioninfo='nii_read_header';
data.fileformat(5).functionread='nii_read_volume';

data.fileformat(6).ext='*.raw';
data.fileformat(6).type='RAW files';
data.fileformat(6).folder='raw';
data.fileformat(6).functioninfo='raw_read_header';
data.fileformat(6).functionread='raw_read_volume';

data.fileformat(7).ext='*.v3d';
data.fileformat(7).type='V3D Philips Scanner';
data.fileformat(7).folder='v3d';
data.fileformat(7).functioninfo='v3d_read_header';
data.fileformat(7).functionread='v3d_read_volume';

data.fileformat(8).ext='*.vmp';
data.fileformat(8).type='VMP BrainVoyager';
data.fileformat(8).folder='vmp';
data.fileformat(8).functioninfo='vmp_read_header';
data.fileformat(8).functionread='vmp_read_volume';

data.fileformat(9).ext='*.xif';
data.fileformat(9).type='XIF HDllab/ATL ultrasound';
data.fileformat(9).folder='xif';
data.fileformat(9).functioninfo='xif_read_header';
data.fileformat(9).functionread='xif_read_volume';

data.fileformat(10).ext='*.vtk';
data.fileformat(10).type='VTK Visualization Toolkit';
data.fileformat(10).folder='vtk';
data.fileformat(10).functioninfo='vtk_read_header';
data.fileformat(10).functionread='vtk_read_volume';

data.fileformat(11).ext='*.mha';
data.fileformat(11).type='Insight Meta-Image';
data.fileformat(11).folder='mha';
data.fileformat(11).functioninfo='mha_read_header';
data.fileformat(11).functionread='mha_read_volume';

data.fileformat(12).ext='*.vff';
data.fileformat(12).type='Micro CT';
data.fileformat(12).folder='vff';
data.fileformat(12).functioninfo='vff_read_header';
data.fileformat(12).functionread='vff_read_volume';

data.fileformat(13).ext='*.par';
data.fileformat(13).type='Philips PAR/REC';
data.fileformat(13).folder='par';
data.fileformat(13).functioninfo='par_read_header';
data.fileformat(13).functionread='par_read_volume';

%---- End supported file formats ----%

% Get path of ReadData3D
functionname='ReadData3D.m';
functiondir=which(functionname);
functiondir=functiondir(1:end-length(functionname));

% Add the file-reader functions also to the matlab path
addpath([functiondir '/subfunctions']);
for i=1:length(data.fileformat), addpath(fullfile(functiondir, data.fileformat(i).folder)); end

% Make popuplist file formats
fileformatcell=cell(1,length(data.fileformat));
for i=1:length(data.fileformat), fileformatcell{i}=[data.fileformat(i).type '   (' data.fileformat(i).ext ')']; end
set(handles.popupmenu_format,'String',fileformatcell);

% Check if last filename is present from a previous time
data.configfile=[functiondir '/lastfile.mat'];
filename='';
fileformatid=1;
if(exist(data.configfile,'file')), load(data.configfile); end
data.handles=handles;
data.lastfilename=[];
data.volume=[];
data.info=[];

% If filename is selected, look if the extention is known
found=0;
if(~isempty(varargin))
    filename=varargin{1}; [pathstr,name,ext]=fileparts(filename);
    for i=1:length(data.fileformat)
        if(strcmp(data.fileformat(i).ext(2:end),ext)), found=1; fileformatid=i; end
    end
end

% Rescale the databack to original units.
if(length(varargin)>1), real=varargin{2}; else real=true; end

data.real=real;
data.filename=filename;
data.fileformatid=fileformatid;

set(handles.checkbox_real,'Value',data.real);
set(handles.edit_filename,'String',data.filename)
set(handles.popupmenu_format,'Value',data.fileformatid);

% Store all data
setMyData(data);

if(found==0)
    % Show Dialog File selection
    uiwait(handles.figure1);
else
    % Load the File directly
    loaddata();
end

% --- Outputs from this function are returned to the command line.
function varargout = ReadData3D_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
if(ishandle(hObject))
    data=getMyData();
else
    data=[];
end
if(~isempty(data))
    varargout{1} = data.volume;
    varargout{2} = data.info;
else
    varargout{1}=[];
    varargout{2}=[];
end
if(ishandle(hObject))
    close(hObject)
end

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

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

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

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

% --- Executes on button press in pushbutton_browse.
function pushbutton_browse_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_browse (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
data=getMyData();
[extlist extlistid]=FileDialogExtentionList(data);
[filename, dirname,filterindex] = uigetfile(extlist, 'Select a dicom file',fileparts(data.filename));
if(filterindex>0)
    if(extlistid(filterindex)~=0)
        data.fileformatid=extlistid(filterindex);
        set( handles.popupmenu_format,'Value',data.fileformatid);
    end
    if(filename==0), return; end
    filename=[dirname filename];
    data.filename=filename;
    setMyData(data);
    set(handles.edit_filename,'String',data.filename)
end

function [extlist extlistid]=FileDialogExtentionList(data)
extlist=cell(length(data.fileformat)+1,2);
extlistid=zeros(length(data.fileformat)+1,1);
ext=data.fileformat(data.fileformatid).ext;
type=data.fileformat(data.fileformatid).type;
extlistid(1)=data.fileformatid;
extlist{1,1}=ext; extlist{1,2}=[type ' (' ext ')'];
j=1;
for i=1:length(data.fileformat);
    if(i~=data.fileformatid)
        j=j+1;
        ext=data.fileformat(i).ext;
        type=data.fileformat(i).type;
        extlistid(j)=i;
        extlist{j,1}=ext; extlist{j,2}=[type ' (' ext ')'];
    end
end
extlist{end,1}='*.*';
extlist{end,2}='All Files (*.*)';

function setMyData(data)
% Store data struct in figure
setappdata(gcf,'dataload3d',data);

function data=getMyData()
% Get data struct stored in figure
data=getappdata(gcf,'dataload3d');

% --- Executes on button press in pushbutton_cancel.
function pushbutton_cancel_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_cancel (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
setMyData([]);
uiresume;

% --- Executes on selection change in popupmenu_format.
function popupmenu_format_Callback(hObject, eventdata, handles)
% hObject    handle to popupmenu_format (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 popupmenu_format contents as cell array
%        contents{get(hObject,'Value')} returns selected item from popupmenu_format
data=getMyData();
data.fileformatid=get( handles.popupmenu_format,'Value');
setMyData(data);

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

% Hint: popupmenu 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 pushbutton_load.
function pushbutton_load_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_load (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
data=getMyData();
data.filename=get(handles.edit_filename,'string');
loaddata();
pause(0.1);
uiresume

function loaddata()
data=getMyData();
set(data.handles.figure1,'Pointer','watch'); drawnow('expose');
if(~strcmp(data.lastfilename,data.filename))
    % Get info
    fhandle = str2func( data.fileformat(data.fileformatid).functioninfo);
    data.info=feval(fhandle,data.filename);
    data.lastfilename=data.filename;
end
fhandle = str2func( data.fileformat(data.fileformatid).functionread);
data.volume=feval(fhandle,data.info);
if(data.real)
    data.volume=single(data.volume);
    if(isfield(data.info,'RescaleSlope')),
        data.volume=data.volume*data.info.RescaleSlope;
    else
        disp('RescaleSlope not available, assuming 1')
    end
    if(isfield(data.info,'RescaleIntercept')),
        data.volume=data.volume+data.info.RescaleIntercept;
    else
        disp('RescaleIntercept not available, assuming 0')
    end
end
setMyData(data);
set(data.handles.figure1,'Pointer','arrow')

% Save the filename, for the next time this function is used
filename=data.filename; fileformatid=data.fileformatid;
try save(data.configfile,'filename','fileformatid'); catch ME; disp(ME.message); end

% --- Executes on button press in pushbutton_info.
function pushbutton_info_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_info (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
data=getMyData();
data.filename=get(handles.edit_filename,'string');
if(~strcmp(data.lastfilename,data.filename))
    % Get info
    set(data.handles.figure1,'Pointer','watch'); drawnow('expose');
    fhandle = str2func( data.fileformat(data.fileformatid).functioninfo);
    data.info=feval(fhandle,data.filename);
    data.lastfilename=data.filename;
    set(data.handles.figure1,'Pointer','arrow')
end
setMyData(data);
% Show info
InfoData3D(data.info);

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

% Hint: get(hObject,'Value') returns toggle state of checkbox_real
data=getMyData();
data.real=get(handles.checkbox_real,'Value');
setMyData(data);

author: tomdoel

repo: pulmonarytoolkit


text:
function [costPatchCand, uvBiasCand] = ...
    sc_patch_cost(trgPatch, srcPatch, modelPlane, uvPlaneIDData, ...
    trgPos, srcTform, srcPosMap, bdPos, optS)

numUvPix      = size(srcPatch, 3);
costPatchCand = zeros(numUvPix, 4, 'single');

% Patch cost - appearance
[costApp, uvBiasCand] = sc_patch_cost_app(trgPatch, srcPatch, optS);

% Patch cost - planar guidance
srcPos = srcTform(:,7:8);
costPlane = sc_patch_cost_plane(modelPlane.mLogLPlaneProb, uvPlaneIDData, trgPos, srcPos);

if(optS.useCoherence) % Patch cost - coherence
    costCoherence = sc_patch_cost_coherence(trgPos, srcPosMap, srcTform, optS);
end

% Patch cost - directional guidance
costDirect = sc_patch_cost_direct(uvPlaneIDData, trgPos, srcPos, modelPlane, optS);

% Patch cost - proximity cost
costProx = sc_patch_cost_prox(srcPos, trgPos, bdPos, optS);

% costReg = double(uvPixUpdateSrc == 1);

% Weighted sum of the appearance cost and guidance cost
costPatchCand(:,1) = costApp;
costPatchCand(:,2) = optS.lambdaPlane * costPlane;

if(optS.useCoherence)
    costPatchCand(:,3) = optS.lambdaCoherence * costCoherence;
end

costPatchCand(:,4) = optS.lambdaDirect * costDirect;
costPatchCand(:,5) = optS.lambdaProx   * costProx;
% costPatchCand(3,:) = optS.lambdaReg    * costReg;

% [To-Do] Adaptive weighting for patch cost
% costPatchCand(2:5,:) = (iLvl/optS.numPyrLvl)*costPatchCand(2:5,:);

end

function [costApp, uvBias] = sc_patch_cost_app(trgPatch, srcPatch, optS)

% SC_PATCH_COST_APP
%
% Compute the weighted sum of the absolute difference between
% cost between source and target patches

% Input:
%   - trgPatch:  target patch
%   - srcPatch:  source patch
%   - optS:      parameter
% Output:
%   - costApp:   appearance cost
%   - uvBias     bias for the closet neighbor

% =========================================================================
% Apply bias correction
% =========================================================================

uvBias = [];
if(optS.useBiasCorrection)
    % Mean of source and target patch
    meanTrgPatch = mean(trgPatch, 1);
    meanSrcPatch = mean(srcPatch, 1);

    % Compute bias and clamp it to inteval [optS.minBias, optS.maxBias]
    uvBias = meanTrgPatch - meanSrcPatch;
    uvBias = sc_clamp(uvBias, optS.minBias, optS.maxBias);

    % Update the UV map for bias
    %     uvBias = reshape(biasPatch, 3, numUvValidPix);

    % Apply the bias correction to source patch
    srcPatch = bsxfun(@plus, srcPatch, uvBias);
end

% =========================================================================
% Compute weighted distance
% =========================================================================
% Compute distance
patchDist = trgPatch - srcPatch;

% Sum of absolute distance
if(strcmp(optS.costType, 'L1'))
    patchDist = abs(patchDist);
elseif(strcmp(optS.costType, 'L2'))
    patchDist = patchDist.^2;
end

% Apply weights
patchDist = bsxfun(@times, patchDist, optS.wPatch);
costApp  = squeeze(sum(sum(patchDist, 1),2));
end

function cost = sc_patch_cost_coherence(trgPos, srcPosMap, srcTform, opt)

% SC_PATCH_COST_COHERENCE: spatial coherence cost
%
% Input
%   - trgPos:    target patch positions [numPix] x [2]
%   - srcPosMap: source patch map [imgH] x [imgW] x [2]
%   - srcPos:    source patch positions [numPix] x [2]
%   - srcTfmG:   source patch geometric transformation [numPix] x [9]
%   - opt:       parameters
% Output:
%   - cost:      spatio-temporal coherence cost

numPix = size(trgPos,1);
[imgH, imgW, ~] = size(srcPosMap);

% initialize coherence cost
cost = zeros(numPix, 1, 'single');

% =========================================================================
% Compute spatial coherence cost
% =========================================================================
for i = 1:4
    % source patch position prediction
    v = opt.propDir(i,:);
    srcPosP  = sc_trans_tform(srcTform, v);
    srcPosP =  srcPosP(:,7:8);

    % source patch positions of neighboring target patches
    trgPosN = bsxfun(@plus, trgPos, v);
    trgIndN = uint32(sub2ind([imgH, imgW], trgPosN(:,2), trgPosN(:,1)));
    srcPosN = sc_uvMat_from_uvMap(srcPosMap,  trgIndN);

    validSrc = (srcPosN(:,1) ~= 0);

    % add cost if the differences are high
    diff           = zeros(numPix, 1, 'single');
    diff(validSrc) = sum(abs(srcPosP(validSrc,:) - srcPosN(validSrc,:)), 2) > 1;
    cost    = cost + diff;
end

% cost = opt.lambdaCoherence*cost;

end

function costPlane = sc_patch_cost_plane(mLogLPlaneProb, uvPlaneIDData, trgPixSub, srcPixSub)

% SC_PATCH_COST_PLANE
%
% Compute planar costs (See Eqn 11 in the paper)

[imgH, imgW, numPlane] = size(mLogLPlaneProb);

srcPixSub     = round(srcPixSub);
uvPlaneIDData = double(uvPlaneIDData);

trgPixIndCur = sub2ind([imgH, imgW, numPlane], trgPixSub(:,2), trgPixSub(:,1), uvPlaneIDData);
srcPixIndCur = sub2ind([imgH, imgW, numPlane], srcPixSub(:,2), srcPixSub(:,1), uvPlaneIDData);

costPlane = mLogLPlaneProb(trgPixIndCur) + mLogLPlaneProb(srcPixIndCur);

end

function costProx = sc_patch_cost_prox(srcPos, trgPos, uvDtBdPixPos, optS)

% SC_PATCH_COST_PROX

% Encourage to sample patches near to the target patch

d = srcPos - trgPos;
d = sqrt(sum(d.^2,2));

d = d/optS.imgSize;
uvDtBdPixPos = uvDtBdPixPos/optS.imgSize;

% costProx = max(0, d - optS.proxThres);
% Shrinkage thresholding
costProx = max(0, d - uvDtBdPixPos - optS.proxThres);
% costProx = max(0, d - optS.proxThres);

end

function costDirect = sc_patch_cost_direct(uvPlaneIDData, trgPos, srcPos, modelPlane, optS)

% SC_PATCH_COST_DIRECT
%
% Compute the directional cost (See Eqn 13 in the paper)
%
% Input
%   -
% Output
%   -

numUvPix = size(trgPos, 1);
costDirect = optS.lambdaDirect*ones(numUvPix, 2, 'single');

for indPlane = 1: modelPlane.numPlane
    % Retrieve the uv pixels that have the current plane label
    uvPlaneIndCur  = uvPlaneIDData == indPlane;
    numPlanePixCur = sum(uvPlaneIndCur);

    if(indPlane == modelPlane.numPlane)
        costDirect(uvPlaneIndCur, :) = optS.imgSize*optS.directThres;
    else
        % The rectifying transformation for the plane
        rectMat = modelPlane.rectMat{indPlane};
        h7 = rectMat(3,1);
        h8 = rectMat(3,2);

        if(numPlanePixCur~=0)
            trgPosCur = trgPos(uvPlaneIndCur, :) - 1;
            srcPosCur = srcPos(uvPlaneIndCur, :) - 1;

            for iTheta = 1:2
                rotMat = modelPlane.rotMat{indPlane, iTheta};

                rotRecMat = rotMat;
                rotRecMat(3,1) = h7;    rotRecMat(3,2) = h8;
                rotRecMat = rotRecMat';

                trgPosCurRect = cat(2, trgPosCur, ones(numPlanePixCur,1))*rotRecMat;
                trgPosCurRect = bsxfun(@rdivide, trgPosCurRect, trgPosCurRect(:,3));

                % Source patch center position in the rectified domain
                srcPosCurRect = cat(2, srcPosCur, ones(numPlanePixCur,1))*rotRecMat;
                srcPosCurRect = bsxfun(@rdivide, srcPosCurRect, srcPosCurRect(:,3));

                costDirect(uvPlaneIndCur, iTheta) = abs(srcPosCurRect(:,2) - trgPosCurRect(:,2));
            end
        end
    end
end

costDirect = min(costDirect, [], 2);
costDirect = min(costDirect/optS.imgSize, optS.directThres);

end

repo: StructCompletion


text:
function track_mha(srcdir, srcpattern, lookahead, load_data, destmat)
%TRACK_MHA Track eddies using Multiple Hypothesis Assignment
% srdir:      Directory in which to find saved eddy frames
% srcpattern: Filename pattern used for eddies (ex. anticyc_ for
%             anticyc_19921014.mat)
% lookahead:  Boolean value to run the scan with lookahead enabled [0 or 1]
% load data:  Path to pre-generated data. Use 0 to load no data.
% destmat:    Path to saved mat-file. (should NOT end in .mat)
    if isa(load_data, 'numeric')
        load_data = '0';
    end

    mha_path = mfilename('fullpath');
    sep = strfind(mha_path, filesep());
    mha_path = [mha_path(1:sep(end)) 'track_mha.py'];

    cmd = sprintf('python "%s" "%s" "%s" %d "%s" "%s"', mha_path, ...
        srcdir, srcpattern, lookahead, load_data, destmat);
    status = system(cmd, '-echo');
    if status
        error('Failed to run MHA');
    end
end

author: jfaghm

repo: OceanEddies

save_path: github-repos/MATLAB/jfaghm-OceanEddies


text:
function gf=gradloglikGaP(x, varargin)
% f=loglikGaP(x, varargin)
% complete log likellihood of the GaP model
% Vxt = varargin{1};      %data
% sigpsf = varargin{2};  %std deviation of the PSF gaussian approx
% alpha = varargin{3}; %parameters of the Gamma prior on the blinking
% beta = varargin{4}; %parameters of the Gamma prior on the blinking
% peval = varargin{5}; %parameters
% x(1:end-2*peval.ncomp) is Hkt

Vxt = varargin{1};      %data
sigpsf = varargin{2};  %std deviation of the PSF gaussian approx
alpha = varargin{3}; %parameters of the Gamma prior on the blinking
beta = varargin{4}; %parameters of the Gamma prior on the blinking
peval = varargin{5}; %parameters

[Hkt, cx, cy, Wxk] = loglikGaPreadparam(x,varargin);

[Wxkbg,Hktbg]=addbg(Wxk, Hkt, peval.bg);
Hkt=max(Hkt,eps); % adjust small values to avoid undeflow
P=Wxkbg*Hktbg; %current approximation

%linear grasdient shifted by cx
xxvc = lineargrad([peval.nx, peval.ny, peval.ncomp], cx, 'xx');
yyvc = lineargrad([peval.nx, peval.ny, peval.ncomp], cy, 'yy');

% dW/dcx:
Wxtcx=1/sigpsf^2*xxvc.*Wxk;
% dW/dcy:
Wxtcy=1/sigpsf^2*yyvc.*Wxk;

% d(log(L))/dHkt:
% gfHkt=(alpha-1)*1./Hkt - 1/beta*Hkt + Wxk'*(Vxt./P)-1;
gfHkt= Wxk'*(Vxt./P)-1; %without background (->not Wxkgb) and d(log(L)/dcx)
% d(log(L))/dcx:
gfcx=diag(Wxtcx'*(Vxt./P-ones(peval.nx*peval.ny, peval.nt))*Hkt');
% d(log(L))/dcy:
gfcy=-diag(Wxtcy'*(Vxt./P-ones(peval.nx*peval.ny, peval.nt))*Hkt');

% gf = [reshape(gfHkt,1,peval.nt*peval.ncomp), gfcx', gfcy'];
gf = [reshape(gfHkt,1,peval.nt*peval.ncomp)];
% gf = [gfcx', gfcy'];
gf=-gf; %conjugate gradient is minimizing!
end
function Wnorm=normalizePSF(W)
sw=size(W);
Wr=reshape(W, sw(1)*sw(2),sw(3));
q=squeeze(sum(Wr,1));
Wrnorm=Wr./repmat(q,sw(1)*sw(2),1);
Wnorm=reshape(Wrnorm,sw(1), sw(2), sw(3));
end
function xxvc = lineargrad(sizevec, cx, dir)
switch dir
    case 'xx'
%         xxp=double(xx(sizevec, 'true')); %linear function - pixels
        xxp=double(xx(sizevec, 'corner')); %linear function - pixels
    case 'yy'
%         xxp=double(yy(sizevec, 'true')); %linear function - pixels
        xxp=double(yy(sizevec, 'corner')); %linear function - pixels
    otherwise
        error('Wrong dir')
end

xxv=reshape(xxp,sizevec(1)*sizevec(2),sizevec(3)); %linear function - vector
xxvc=xxv-repmat(cx,sizevec(1)*sizevec(2),1);
% xxp=double(xx([peval.nx, peval.ny, peval.ncomp], 'true')); %linear function - pixels
%yyp=double(yy([peval.nx, peval.ny, peval.ncomp], 'true'));
% xxv=reshape(xxp,peval.nx*peval.ny,peval.ncomp); %linear function - vector
%yyv=reshape(yyp,peval.nx*peval.ny,peval.ncomp);
% xxvc=xxv-repmat(cx,peval.ncomp,1);
%yyvc=yyv-repmat(cy,peval.ncomp,1);
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
 function c = horzcat(varargin)
%function c = horzcat(varargin)
% called for c = [a, b, ...] where any of them is a fatrix

c = block_fatrix(varargin, 'type', 'row');

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
#!/usr/bin/octave
% Copyright 2018-2020 Marc René Schädler
%
% This file is part of the mobile hearing aid prototype project
% The the mobile hearing aid prototype project 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.
%
% The mobile hearing aid prototype project 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 the mobile hearing aid prototype project. If not, see http://www.gnu.org/licenses/.

cd(fileparts(mfilename('fullpath')))

close all
clear
clc

graphics_toolkit qt;

thresholds_freqs = [125 250 500 1000 2000 4000 8000 16000];
thresholds_init  = [ 39  71  71   62   60   58   58    36]; % dB SPL at mic
normal_threshold = [ 33  24  13    5   13   15   14    11]; % Estimated

plot_freqs = 1000 .* 2.^(-3:0.5:4);
plot_threshold = interp1(thresholds_freqs,normal_threshold,plot_freqs);
plot_levels = 20:15:110;
plot_colors = lines(length(plot_levels));

set (0, "defaultaxesfontsize", 20);
set (0, "defaulttextfontsize", 20);
h_figure = figure('Position', [ 500 200 1500 1000],
       'NumberTitle', 'off',
       'Name', 'Simple self-fitting GUI',
       'toolbar', 'none',
       'menubar', 'none');
h.thresholds_freqs = thresholds_freqs;
h.plot_freqs = plot_freqs;
h.plot_threshold = plot_threshold;
h.plot_levels = plot_levels;
h.plot_colors = plot_colors;

figure_positions = {[0.075 0.37 0.4 0.5],[0.55 0.37 0.4 0.5]};
position_title = {'L E F T', 'R I G H T'};
for i=1:2
  h.ax(i) = axes ('position', figure_positions{i});
  box('on');
  h_threshold = plot(log(thresholds_freqs), zeros(size(thresholds_freqs)),'-k','linewidth',3);
  hold on;
  title(position_title{i});
  plot(log(plot_freqs), plot_threshold,'--k','linewidth',3);
  xlim(log([min(plot_freqs) max(plot_freqs)]));
  ylim([0 130]);
  set(gca,'xtick',log(plot_freqs(1:2:end)));
  set(gca,'xticklabel',plot_freqs(1:2:end));
  set(gca,'ytick',0:10:100);
  xlabel('Frequency / Hz');
  ylabel('Levels / dB SPL (in device)');
  for j=1:length(plot_levels)
    plot(log(plot_freqs),ones(size(plot_freqs)).*plot_levels(j),'--','color',plot_colors(j,:),'linewidth',3);
  end
  h_amplified = zeros(size(plot_levels));
  for j=1:length(plot_levels)
    h_amplified(j) = plot(log(plot_freqs),ones(size(plot_freqs)).*plot_levels(j),'-','color',plot_colors(j,:),'linewidth',3);
  end
  h.h_activedata{i} = h_amplified;
  h.h_threshold{i} = h_threshold;
end

function s = gt_data2string(gt_data)
  [num_rows, num_cols] = size(gt_data);
  s = '[';
  for i=1:num_rows
    s = [s, '['];
    for j=1:num_cols
      s = [s, sprintf('%.2f',gt_data(i,j))];
      if j < num_cols
        s = [s, ' '];
      else
        s = [s, ']'];
      end
    end
    if i < num_rows
      s = [s, ';'];
    else
      s = [s, ']'];
    end
  end
end

function text2speech(message)
  system(['echo text2speech "',message,'" | nc -w 1 127.0.0.1 33338']);
end

function mhacontrol(command)
  system(['echo mhacontrol "',command,'" | nc -w 1 127.0.0.1 33338']);
end

function mhaplay(filename)
  system(['echo mhaplay "',filename,'" "yes" | nc -w 1 127.0.0.1 33338']);
end

function thresholdnoise(status)
  system(['echo thresholdnoise "',status,'" | nc -w 1 127.0.0.1 33338']);
end

function live(status)
  system(['echo live "',status,'" | nc -w 1 127.0.0.1 33338']);
end

function feedback(duration)
  if isnumeric(duration)
    duration = num2str(duration);
  end
  system(['echo feedback "',duration,'" | nc -w 1 127.0.0.1 33338']);
end

function record(duration)
  system(['echo record "',duration,'" | nc -w 1 127.0.0.1 33338']);
end

function send_gaintable(gt_data)
  system(['echo " | nc -w 1 127.0.0.1 33338']);
end

function update_gaintable (obj)
  tic
  h = guidata (obj);
  offset = 30+30.*(0.5-get(h.offset_slider,'value'));
  rolloff = 1+2.^-(get(h.rolloff_slider,'value').*3);
  marginfactor = get(h.marginfactor_slider,'value');
  center = 50+get(h.center_slider,'value').*40;
  focus = round(get(h.focus_slider,'value').*100)./10;

  thresholds_left = zeros(size(h.thresholds_freqs));
  thresholds_right = zeros(size(h.thresholds_freqs));

  for i=1:length(h.thresholds_freqs)
    thresholds_left(i) = str2num(get(h.thresholds_left{i},'string'));
    thresholds_right(i) = str2num(get(h.thresholds_right{i},'string'));
  end

  set(h.h_threshold{1},'ydata',thresholds_left);
  set(h.h_threshold{2},'ydata',thresholds_right);

  [gt_data, gt_freqs, gt_levels] = prescription_minimalistic(h.thresholds_freqs, thresholds_left, thresholds_right, offset, marginfactor, rolloff, center, focus);
  for i=1:2
    gain = interp2(gt_levels,gt_freqs.',gt_data(1+(i-1).*length(gt_freqs):i.*length(gt_freqs),:),h.plot_levels.',h.plot_freqs,'linear');
    for j=1:length(h.plot_levels)
      set(h.h_activedata{i}(j),'ydata',h.plot_levels(j)+gain(:,j));
    end
  end
  mhacontrol(['mha.transducers.mhachain.altplugs.dynamiccompression.mhachain.dc.gtdata = ',gt_data2string(gt_data)]);
  toc
  drawnow;
end

function update_live(obj)
  h = guidata (obj);
  if get(h.live_checkbox,'value')
    live('on');
  else
    live('off');
  end
end

function update_noise(obj)
  h = guidata (obj);
  if get(h.noise_checkbox,'value')
    thresholdnoise('on');
  else
    thresholdnoise('off');
  end
end

function start_recording()
  record('10');
end

function update_amplification(obj)
  h = guidata (obj);
  if get(h.amplification_checkbox,'value')
    mhacontrol('mha.transducers.mhachain.altplugs.select = dynamiccompression');
  else
    mhacontrol('mha.transducers.mhachain.altplugs.select = identity');
  end
end

function update_calib(obj)
  h = guidata (obj);
  if get(h.calib_checkbox,'value')
	% values taken from current openMHA.cfg
    mhacontrol('mha.transducers.calib_in.fir = [[-0.674757033017875 -0.339128228007168 0.030117856798366 -0.119542483672035 0.131726684555067 -0.071387589553860 0.056561094657502 -0.013365060531109 0.007091265106959 0.006913350591302 -0.002602409570910 0.023338189483743 -0.018652297865493 0.021281416315377 -0.014374905102702 0.015659551668403];[-0.674757033017875 -0.339128228007168 0.030117856798366 -0.119542483672035 0.131726684555067 -0.071387589553860 0.056561094657502 -0.013365060531109 0.007091265106959 0.006913350591302 -0.002602409570910 0.023338189483743 -0.018652297865493 0.021281416315377 -0.014374905102702 0.015659551668403]]');
    mhacontrol('mha.transducers.calib_out.fir = [[-0.723329250680037 1.285750531421244 -1.764370643350338 1.702299222653803 0.569187694811727 -1.301642333817312 0.540112043557634 -0.111972475691304 0.187607100395961 0.532899301077540 -0.733701264924906 0.290587502735630 0.470074681167595 -0.330332208101367 -0.010615612830630 -0.221018479662762 -0.020491963942295 0.385847945936633 -0.198969798086483 -0.123437157204177 0.152968427130900 -0.205009361694624 -0.071974598348291 0.041537565470939 -0.058703551389704 0.006171609199330 -0.037439953306526 -0.032735557865821 0.028903184718236 -0.019431703287974 -0.069837180982837 0.038967963643885 -0.015801874625706 -0.060091628384134 0.057102986709234 -0.043134256115441 -0.016815289718170 0.002879163232524 -0.031134102833679 0.026124175604008 -0.040758325083950 0.022598560658844 0.012982112833010 -0.029909215215110 0.018353193685970 0.012015086827306 0.031537385134343 -0.028980188877489 0.026113358679986 0.017332696993042 -0.018853690496196 0.059665425003585 -0.017281793549857 0.007285675861939 0.010218784848149 0.003173184622963 0.040817408121183 -0.024810702312558 0.007344195504224 0.007639434860116 -0.004797394889051 0.025321400227842 -0.005638659121956 -0.002191206231116];[-0.723329250680037 1.285750531421244 -1.764370643350338 1.702299222653803 0.569187694811727 -1.301642333817312 0.540112043557634 -0.111972475691304 0.187607100395961 0.532899301077540 -0.733701264924906 0.290587502735630 0.470074681167595 -0.330332208101367 -0.010615612830630 -0.221018479662762 -0.020491963942295 0.385847945936633 -0.198969798086483 -0.123437157204177 0.152968427130900 -0.205009361694624 -0.071974598348291 0.041537565470939 -0.058703551389704 0.006171609199330 -0.037439953306526 -0.032735557865821 0.028903184718236 -0.019431703287974 -0.069837180982837 0.038967963643885 -0.015801874625706 -0.060091628384134 0.057102986709234 -0.043134256115441 -0.016815289718170 0.002879163232524 -0.031134102833679 0.026124175604008 -0.040758325083950 0.022598560658844 0.012982112833010 -0.029909215215110 0.018353193685970 0.012015086827306 0.031537385134343 -0.028980188877489 0.026113358679986 0.017332696993042 -0.018853690496196 0.059665425003585 -0.017281793549857 0.007285675861939 0.010218784848149 0.003173184622963 0.040817408121183 -0.024810702312558 0.007344195504224 0.007639434860116 -0.004797394889051 0.025321400227842 -0.005638659121956 -0.002191206231116]]');
  else
    mhacontrol('mha.transducers.calib_in.fir = 1');
    mhacontrol('mha.transducers.calib_out.fir = 1');
  end
end

function measure_feedback(obj)
  feedback(3);
end

function update_playback(obj)
  h = guidata (obj);
  string = get(h.playback_popup,'string');
  value = get(h.playback_popup,'value');
  selection = string{value};
  switch selection
    case 'off'
      filename = '';
    case 'last record'
      filename ='/dev/shm/recording.wav';
    otherwise
      filename = ['/home/pi/hearingaid-prototype/recordings/' selection '.wav'];
  end
  mhaplay(filename);
end

for i=1:length(thresholds_freqs)
  uicontrol ('style', 'text',
    'units', 'normalized',
    'string', num2str(thresholds_freqs(i)) ,
    'horizontalalignment', 'center',
    'position', [0.02+i*0.05 0.95 0.05 0.05]);
  h.thresholds_left{i} = uicontrol ('style', 'edit',
    'units', 'normalized',
    'string', num2str(thresholds_init(i)),
    'callback', @update_gaintable,
    'position', [0.02+i*0.05 0.90 0.05 0.05]);
end

for i=1:length(thresholds_freqs)
  uicontrol ('style', 'text',
    'units', 'normalized',
    'string', num2str(thresholds_freqs(i)) ,
    'horizontalalignment', 'center',
    'position', [0.50+i*0.05 0.95 0.05 0.05]);
  h.thresholds_right{i} = uicontrol ('style', 'edit',
    'units', 'normalized',
    'string', num2str(thresholds_init(i)),
    'callback', @update_gaintable,
    'position', [0.50+i*0.05 0.90 0.05 0.05]);
end

uicontrol ('style', 'text',
  'units', 'normalized',
  'string', 'gain' ,
  'horizontalalignment', 'left',
  'position', [0.05 0.25 0.15 0.05]);

h.offset_slider = uicontrol ('style', 'slider',
  'units', 'normalized',
  'sliderstep', [0.1 0.1],
  'string', 'slider',
  'callback', @update_gaintable,
  'value', 0.5,
  'position', [0.15 0.25 0.65 0.05]);

uicontrol ('style', 'text',
  'units', 'normalized',
  'string', 'compression' ,
  'horizontalalignment', 'left',
  'position', [0.05 0.2 0.15 0.05]);

h.rolloff_slider = uicontrol ('style', 'slider',
  'units', 'normalized',
  'sliderstep', [0.1 0.1],
  'string', 'slider',
  'callback', @update_gaintable,
  'value', 0.5,
  'position', [0.15 0.2 0.65 0.05]);

uicontrol ('style', 'text',
  'units', 'normalized',
  'string', 'boost factor' ,
  'horizontalalignment', 'left',
  'position', [0.05 0.15 0.15 0.05]);

h.marginfactor_slider = uicontrol ('style', 'slider',
  'units', 'normalized',
  'sliderstep', [0.1 0.1],
  'string', 'slider',
  'callback', @update_gaintable,
  'value', 0.5,
  'position', [0.15 0.15 0.65 0.05]);

uicontrol ('style', 'text',
  'units', 'normalized',
  'string', 'boost level' ,
  'horizontalalignment', 'left',
  'position', [0.05 0.1 0.15 0.05]);

h.center_slider = uicontrol ('style', 'slider',
  'units', 'normalized',
  'sliderstep', [0.1 0.1],
  'string', 'slider',
  'callback', @update_gaintable,
  'value', 0.5,
  'position', [0.15 0.1 0.65 0.05]);

uicontrol ('style', 'text',
  'units', 'normalized',
  'string', 'focus' ,
  'horizontalalignment', 'left',
  'position', [0.05 0.05 0.15 0.05]);

h.focus_slider = uicontrol ('style', 'slider',
  'units', 'normalized',
  'sliderstep', [0.1 0.1],
  'string', 'slider',
  'callback', @update_gaintable,
  'value', 0.5,
  'position', [0.15 0.05 0.65 0.05]);

h.live_checkbox = uicontrol ('style', 'checkbox',
  'units', 'normalized',
  'string', 'live',
  'value', 1,
  'callback', @update_live,
  'position', [0.85 0.275 0.05 0.025]);

h.record_button = uicontrol ('style', 'pushbutton',
  'units', 'normalized',
  'string', 'record',
  'callback', @start_recording,
  'position', [0.90 0.275 0.05 0.025]);

uicontrol ('style', 'text',
  'units', 'normalized',
  'string', 'playback' ,
  'horizontalalignment', 'left',
  'position', [0.85 0.25 0.1 0.025]);

h.playback_popup = uicontrol ('style', 'popupmenu',
  'units', 'normalized',
  'string', {'off' 'last record' '35-29-25_fridge' '36-26-19_clock' '53-41-34_footsteps' '55-44-33_handwash' '63-57-37_microwave' '63-57-51_coffeemachine' '64-33-24_keyboardwriting' '64-42-35_foyer' '64-53-46_autobahn' '66-56-48_siren' '68-57-50_drivingnews' '68-59-52_street' '70-61-53_street' '72-57-49_bird' '72-63-60_coffeegrinder' '73-50-42_steeldrum-steps' '73-60-50_cardriveby' '73-63-57_steet-laughter' '75-25-17_phone' '77-54-46_ducks' '78-68-60_street' '78-69-62_lisboa-bride-bird' '79-59-50_bar' '81-54-34_coffeemachine' '81-61-47_traincrossing' '81-66-60_lisboa-plane-landing' '82-53-35_shakeout' '83-69-60_mainstreet' '83-73-65_traffic' '86-75-68_lisboa-bride-train' '88-73-65_lisboa-train' '100-64-37_flute'},
  'callback', @update_playback,
  'position', [0.85 0.225 0.1 0.025]);

h.amplification_checkbox = uicontrol ('style', 'checkbox',
  'units', 'normalized',
  'string', 'amplification',
  'value', 0,
  'callback', @update_amplification,
  'position', [0.85 0.15 0.1 0.025]);

h.noise_checkbox = uicontrol ('style', 'checkbox',
  'units', 'normalized',
  'string', 'threshold noise',
  'value', 0,
  'callback', @update_noise,
  'position', [0.85 0.175 0.1 0.025]);

h.calib_checkbox = uicontrol ('style', 'checkbox',
  'units', 'normalized',
  'string', 'calibration',
  'value', 1,
  'callback', @update_calib,
  'position', [0.85 0.125 0.1 0.025]);

h.feedback_button = uicontrol ('style', 'pushbutton',
  'units', 'normalized',
  'string', 'feedback',
  'callback', @measure_feedback,
  'position', [0.85 0.05 0.1 0.05]);

guidata(gcf, h);
update_calib(gcf);
update_amplification(gcf);
update_gaintable(gcf);
update_noise(gcf);
update_playback(gcf);
update_live(gcf);
waitfor(h_figure);

author: m-r-s

repo: hearingaid-prototype


text:
function opts = getTrainingOptions(hyperparam,imds_validation)

if isempty(imds_validation) % if there is no validation set
opts = trainingOptions('adam',...
    'InitialLearnRate',hyperparam.InitialLearnRate,...
    'L2Regularization',hyperparam.L2Regularization,...
    'MiniBatchSize',hyperparam.MiniBatchSize,...
    'MaxEpochs',hyperparam.MaxEpochs,...
    'ExecutionEnvironment',hyperparam.ExecutionEnvironment,...
    'VerboseFrequency',50,...
    'Plots','training-progress');
else  % if there is a validation set
opts = trainingOptions('adam',...
    'InitialLearnRate',hyperparam.InitialLearnRate,...
    'ValidationData',imds_validation,...
    'ValidationFrequency',hyperparam.ValidationFrequency,...
    'ValidationPatience',hyperparam.ValidationPatience,...
    'L2Regularization',hyperparam.L2Regularization,...
    'MiniBatchSize',hyperparam.MiniBatchSize,...
    'MaxEpochs',hyperparam.MaxEpochs,...
    'ExecutionEnvironment',hyperparam.ExecutionEnvironment,...
    'VerboseFrequency',50,...
    'Plots','training-progress');
end

end

author: jnkather

repo: MSIfromHE

save_path: github-repos/MATLAB/jnkather-MSIfromHE


%
% The main script.
%
clear all;
addpath('bss_eval');

% sample rate
fs=16000;
% reverberation time
RT60_SET=[0.6, 0.3, 0.8];
% SIR (dB) of the test data
SIR_SET=[0];
% SER (dB) of the test data
SER_SET=[0, -10];
% no. of experiments
numrpts=20;
% for algorithm convergence
UTT_REPEAT=2;
% seed
rngseed=0;

for i=1:length(RT60_SET)
for j=1:length(SIR_SET)
for k=1:length(SER_SET)
    rt60=RT60_SET(i);
    sirdb=SIR_SET(j);
    serdb=SER_SET(k);
    savedir=['output/rt' num2str(rt60) '_sir' num2str(sirdb) '_ser' num2str(serdb)];
    if ~exist(savedir)
        mkdir(savedir);
    end

for rptcount=1:numrpts
    fprintf(1, 'rptcount = %d ...\n', rptcount);

    %% generate test data
    fprintf(1, 'generate test data\n');

    rngseed=(i-1)*1e4+(j-1)*1e3+(k-1)*1e2+rptcount;
    [testdata, target_interf_echo] = simulateDataImageRIR(sirdb, serdb, rt60, UTT_REPEAT, rngseed);
    siglen=size(testdata, 1)/UTT_REPEAT;
    audiowrite([savedir '/input_data_', num2str(rptcount), '.wav'], testdata(end-siglen+1:end, :), fs);
    audiowrite([savedir '/target_interf_echo', num2str(rptcount), '.wav'], target_interf_echo(end-siglen+1:end, :), fs);

    %% perform separation
    nummics=2;
    numrefs=1;

    fprintf(1, '1. WPE-NLMS-BSS\n');
    out=modeWPE_NLMS_BSS(nummics, numrefs, testdata);
    out(size(testdata, 1)+1:end, :)=[];
    audiowrite([savedir '/01_out_wpe_nlms_bss_', num2str(rptcount), '.wav'], out(end-siglen+1:end, :), fs);

    fprintf(1, '2. NLMS-WPE-BSS\n');
    out=modeNLMS_WPE_BSS(nummics, numrefs, testdata);
    out(size(testdata, 1)+1:end, :)=[];
    audiowrite([savedir '/02_out_nlms_wpe_bss_', num2str(rptcount), '.wav'], out(end-siglen+1:end, :), fs);

    fprintf(1, '3. DR-AEC-BSS\n');
    out=modeDR_AEC_BSS(nummics, numrefs, testdata);
    out(size(testdata, 1)+1:end, :)=[];
    audiowrite([savedir '/03_out_dr_aec_bss_', num2str(rptcount), '.wav'], out(end-siglen+1:end, :), fs);

    fprintf(1, '4. AEC-DR-BSS\n');
    out=modeAEC_DR_BSS(nummics, numrefs, testdata);
    out(size(testdata, 1)+1:end, :)=[];
    audiowrite([savedir '/04_out_aec_dr_bss_', num2str(rptcount), '.wav'], out(end-siglen+1:end, :), fs);

    fprintf(1, '5. DRAEC-BSS\n');
    out=modeDRAEC_BSS(nummics, numrefs, testdata);
    out(size(testdata, 1)+1:end, :)=[];
    audiowrite([savedir '/05_out_draec_bss_', num2str(rptcount), '.wav'], out(end-siglen+1:end, :), fs);

    fprintf(1, '6. Joint-SS\n');
    out=modeBSS(nummics, numrefs, testdata);
    out(size(testdata, 1)+1:end, :)=[];
    audiowrite([savedir '/06_out_joint_ss_', num2str(rptcount), '.wav'], out(end-siglen+1:end, :), fs);
end

Main_score;

end
end
end


text:
function i4_is_prime_test ( )

%*****************************************************************************80
%
%% I4_IS_PRIME_TEST tests I4_IS_PRIME.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4_IS_PRIME_TEST\n' );
  fprintf ( 1, '  I4_IS_PRIME reports whether an integer is prime.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  I     I4_IS_PRIME(I)\n' );
  fprintf ( 1, '\n' );

  for i = -2 : 25
    fprintf ( 1, '  %6d  %1d\n', i, i4_is_prime ( i ) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m