text:
function [Population,Groups,K] = spiltVariables(Problem,BU,BD,s)
% Decision variable grouping with RDG2 in SACCEAMII

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

    %% Initialize
    D = length(BU);
    Groups     = zeros(1,D);
    Population = [];
    Separable  = [];
    NonSeparable = [];

    %% Recursive Decomposition Grouping Methods
    xLL  = BD;
    xNew = Problem.Evaluation(xLL);
    yLL  = xNew.objs;
    Population = [Population,xNew];
    X1   = 1;
    X2   = 2:D;
    % Process
    while ~isempty(X2)
        % Interaction detection
        [sub1,Population] = INTERACT(Problem,Population,X1,X2,xLL,yLL,BU,BD);
        if numel(sub1) == numel(X1) % isequal(sub1,X1)
            if numel(X1) == 1
                Separable = [Separable,X1];
            else
                NonSeparable = [NonSeparable,X1];
            end
			X1 = X2(1);
			X2 = X2(2:end);
        else
            X1 = sub1;
            X2 = X2(~ismember(X2,X1));
        end

        if isempty(X2)
            if numel(X1) > 1
                NonSeparable = [NonSeparable,X1];
            else
                Separable = [Separable,X1];
            end
        end
    end

    %% Split the separable decision varables
     if ~isempty(Separable)
         num = floor(numel(Separable)/s);
         if num <= 1
             K   = 1;
             Groups(Separable) = K;
         else
             K = num;
             idx = 1:1:numel(Separable);
             for i = 1 : num-1
                 select = randperm(numel(idx),s);
                 Groups(Separable(idx(select))) = i;
                 idx(select) = [];
             end
             Groups(Separable(idx)) = num;
         end
         if ~isempty(NonSeparable)
             K = K + 1;
             Groups(NonSeparable) = K;
         end
     else
         K = 1;
         Groups(NonSeparable) = K;
     end
end

function [X1,Population] = INTERACT(Problem,Population,X1,X2,xLL,yLL,BU,BD)
    % Determine gamma
    muM   = eps/2;
    gamma = @(n)(n.*muM)./(1-n.*muM);

    xUL     = xLL;
    xUL(X1) = BU(X1);

    % Calculate delta1
    tNew       = Problem.Evaluation(xUL);
    delta1     = yLL - tNew.objs;
    Population = [Population,tNew];

    % Calculate delta2
    xLM     = xLL;
    xLM(X2) = (BU(X2)+BD(X2))/2;
    xUM     = xUL;
    xUM(X2) = (BU(X2)+BD(X2))/2;
    tNew1   = Problem.Evaluation(xLM);
    tNew2   = Problem.Evaluation(xUM);
    delta2  = tNew1.objs - tNew2.objs;
    Population = [Population,tNew1,tNew2];
    % Update
    F4 = [yLL,tNew.objs,tNew1.objs,tNew2.objs];
    threshold = gamma(numel(BU).^0.5+2)*sum(abs(F4));
    if abs(delta1-delta2) > threshold
        if numel(X2) == 1
            X1 = union(X1,X2);
        else
            % Divide X2 into equally-sized groups G1, G2
            mid = floor(length(X2)/2);
            G1  = X2(1:mid);
            G2  = X2(mid+1:end);
            [subX1,Population] = INTERACT(Problem,Population,X1,G1,xLL,yLL,BU,BD);
            [subX2,Population] = INTERACT(Problem,Population,X1,G2,xLL,yLL,BU,BD);
            X1  = union(subX1,subX2);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
filename = 'CantileverSquareSmall';%'CantileverSquareNewFine';%'LshapeTriFine';%'CantileverSquareNew';'ArchTriFine';%'CantileverSquare';%'ArchTriFine';%'CantileverSquare';%'Lshape';%'LshapeTriFine';%'CantileverSquareSmall';%'';%'Lshape';'CantileverSquare';'Lshape';%'CantileverSquare';%'LshapeFine';%'Bridge_quad_coarse';%'Arch_quad_coarse';%'BridgeCool_Quadrilateral_Bilinear_Structured_Coarse';%'Bridge';%'CantileverSquareSmall';
ptype = 'MACRO';
initial_case = 'given';
m1 = 0.0101;
m2 = 0.0101;
%cost = {'compliance'};
cost = {'stressNorm'};
%cost = {'stressNorm','compliance'};
%weights = [0.55,0.45];
weights = 1;
constraint = {'volumeConstraint'};
filterType = 'P1';
constraint_case = 'EQUALITY';

Vfrac_initial = 0.3;
optimality_initial = 1e-5;
constr_initial = 1e-5;

Vfrac_final = 0.3;
optimality_final = 1e-5;
constr_final = 1e-5;

stressNormExponent_initial = 2;
stressNormExponent_final = 64;

optimizer = 'DualNestedInPrimal';
optimizerUnconstrained = 'PROJECTED GRADIENT';

designVariable = 'MicroParams';
ub = 0.989;
lb = 0.011;
homegenizedVariablesComputer = 'ByVademecum';
%
%vademecumFileName = 'SuperEllipseQMax';
%vademecumFileName = 'SuperEllipseQ2';
vademecumFileName = 'SuperEllipseQOptAnalytic';
%
% designVariable = 'Density';
% homegenizedVariablesComputer = 'ByInterpolation';
% method = 'SIMPALL';
% materialType = 'ISOTROPIC';

kfrac = 2;
nsteps = 250;

plotting = false;
printing = true;
monitoring = true;
monitoring_interval = 1;
maxiter = 5000;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function r8mat_print ( m, n, a, title )

%*****************************************************************************80
%
%% R8MAT_PRINT prints an R8MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows in A.
%
%    Input, integer N, the number of columns in A.
%
%    Input, real A(M,N), the matrix.
%
%    Input, string TITLE, a title to be printed.
%
  r8mat_print_some ( m, n, a, 1, 1, m, n, title );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ vX, mX ] = SolveLsL1ComplexIrls( mA, vB, lambdaFctr, numIterations )
% ----------------------------------------------------------------------------------------------- %
%[ vX, mX ] = SolveLsL1ComplexIrls( mA, vB, lambdaFctr, numIterations )
% Solves the 0.5 * || A x - b ||_2 + \lambda || x ||_1 problem using Fixed
% Point Iteration Method. The model allows A, b and x to be Complex.
% Input:
%   - mA                -   Model Matrix.
%                           The model matrix.
%                           Structure: Matrix (m X n).
%                           Type: 'Single' / 'Double' (Complex).
%                           Range: (-inf, inf).
%   - vB                -   Input Vector.
%                           The model known data.
%                           Structure: Vector (m X 1).
%                           Type: 'Single' / 'Double' (Complex).
%                           Range: (-inf, inf).
%   - paramLambda       -   Parameter Lambda.
%                           The L1 Regularization parameter.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: (0, inf).
%   - numIterations     -   Number of Iterations.
%                           Number of iterations of the algorithm.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range {1, 2, ...}.
% Output:
%   - vX                -   Output Vector.
%                           Structure: Vector (n X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References
%   1.  Wikipedia Fixed Point Iteration Method - https://en.wikipedia.org/wiki/Fixed-point_iteration.
% Remarks:
%   1.  A
% Known Issues:
%   1.  A
% TODO:
%   1.  B
% Release Notes:
%   -   1.0.000     07/11/2016
%       *   First realease version.
% ----------------------------------------------------------------------------------------------- %

paramEps = 1e-9;

% vX  = mAA \ vAb;
vX  = pinv(mA) * vB; %<! Dealing with "Fat Matrix"

mAA = mA' * mA;
vAb = mA' * vB;
% mI = eye(size(mA, 1));

mX = zeros([size(vX, 1), numIterations]);
mX(:, 1) = vX;

for ii = 2:numIterations
    vX = (mAA + (lambdaFctr * diag(1./ (abs(vX) + paramEps)))) \ vAb;

    % Doesn't work (Bad conditioning)
    % mX = diag(abs(vX));
    % vX = (mX * mA') * (((mA * mX * mA.') + (lambdaFctr * mI)) \ vB);

    mX(:, ii) = vX;
end

end

author: RoyiAvital

repo: StackExchangeCodes


text:
%%%  AttitudeControl_Sim
clear
path('./icon/',path);
%Run the model initialization file icon/init.m
Init;

%Constant value
RAD2DEG = 57.2957795;
DEG2RAD = 0.0174533;
%throttle when UAV is hovering
THR_HOVER = 0.609;

%% Initial condition
ModelInit_PosE = [0, 0, 0];
ModelInit_VelB = [0, 0, 0];
ModelInit_AngEuler = [0, 0, 0];
ModelInit_RateB = [0, 0, 0];
ModelInit_Rads = 0;
%% control parameter
%Attitude PID parameters
Kp_YAW_AngleRate = 0.5;
Ki_YAW_AngleRate = 0.01;
Kd_YAW_AngleRate = 0.00;
%integral saturation
Saturation_I_RP_Max = 0.3;
Saturation_I_RP_Min = -0.3;
Saturation_I_Y_Max = 0.2;
Saturation_I_Y_Min = -0.2;
%max control angle,default 35deg
MAX_CONTROL_ANGLE_ROLL = 35;
MAX_CONTROL_ANGLE_PITCH  = 35;
%max control angle rate,rad/s
MAX_CONTROL_ANGLE_RATE_PITCH = 220;
MAX_CONTROL_ANGLE_RATE_ROLL = 220;
MAX_CONTROL_ANGLE_RATE_Y = 200;
%% run simulink model
AttitudeControl_Sim

author: RflySim

repo: RflyExpCode

save_path: github-repos/MATLAB/RflySim-RflyExpCode


text:
function test_interp_2d_test03 ( )

%*****************************************************************************80
%
%% TEST_INTERP_2D_TEST03 plots each grid.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 July 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST_INTERP_2D_TEST03\n' );
  fprintf ( 1, '  Plot each grid.\n' );

  g_num = g00_num ( );
  filename = 'grid_00.png';

  for gi = 1 : g_num

    gt = g00_title ( gi );
    gn = g00_size ( gi );

    [ gx, gy ] = g00_xy ( gi, gn );

    plot ( gx, gy, 'b.', 'Markersize', 25 );
    title ( gt );
    grid on
    xlabel ( '<--- X --->' );
    ylabel ( '<--- Y --->' );
    axis ( 'equal' );

    filename = filename_inc ( filename );
    print ( '-dpng', filename );
    fprintf ( 1, '  Created file "%s".\n', filename );

  end

  return
end


text:
% load all result variables

width = 24;
height = 12;
fontsize = 9;
plot_filename = 'pie_household_02';
plot_folder = 'projects/seminar/';

experiments = {'app_dishwasher_r_0.1', ...
                'app_freezer_r_0.2', ...
                'app_fridge_r_0.1', ...
                'app_kettle_r_0.1', ...
                'app_stereo_r_0.2', ...
                'app_stove_r_1', ...
                'app_tv_r_0.2', ...
            };

inferred = zeros(1, length(experiments) + 4);
actual = zeros(1, length(experiments) + 3);
appliance_names = {};
for i = 1:length(experiments)
    experiment = experiments{i};

    filename = ['results/details/weiss/weiss_initial/', experiment, '/default/result1.mat'];
    load(filename);
    setup = result.setup;
    household = setup.household;
    dataset = setup.dataset;
    granularity = setup.granularity;

    appliance_name = result.appliance_names{1};
    appliance_names = [appliance_names, appliance_name];
    % inferred consumption
    inferred(i) = sum(result.consumption) / 3600;

    evaluation_days = result.evaluation_and_training_days{1};
    appliance_id = getApplianceID(appliance_name);
    appliance_consumption = read_plug_data(dataset, ...
                                    household, ...
                                    appliance_id, ...
                                    evaluation_days, ...
                                    granularity);
    actual_valid = appliance_consumption ~= -1;
    percentage_valid = sum(actual_valid) / length(appliance_consumption);
    actual(i) = sum(appliance_consumption(actual_valid)) / 3600 / percentage_valid;
end

total_consumption = read_smartmeter_data(dataset, ...
                                    num2str(household, '%02d'), ...
                                    evaluation_days, ...
                                    granularity, ...
                                    'powerallphases');

%% Laptop
num_days = size(evaluation_days,1);
laptop = num_days * 40 * 4;
idx = length(experiments) + 1;
inferred(idx) = laptop;
appliance_id = getApplianceID('Laptop');
appliance_consumption = read_plug_data(dataset, ...
                                household, ...
                                appliance_id, ...
                                evaluation_days, ...
                                granularity);
actual_valid = appliance_consumption ~= -1;
percentage_valid = sum(actual_valid) / length(appliance_consumption);
actual(idx) = sum(appliance_consumption(actual_valid)) / 3600 / percentage_valid;
appliance_names{idx} = 'Laptop';

%% Lamp
lamp = num_days * 100 * 1;
idx = length(experiments) + 2;
inferred(idx) = lamp;
appliance_id = getApplianceID('Lamp');
appliance_consumption = read_plug_data(dataset, ...
                                household, ...
                                appliance_id, ...
                                evaluation_days, ...
                                granularity);
actual_valid = appliance_consumption ~= -1;
percentage_valid = sum(actual_valid) / length(appliance_consumption);
actual(idx) = sum(appliance_consumption(actual_valid)) / 3600 / percentage_valid;
appliance_names{idx} = 'Lamp';

%% Standby
num_days = size(evaluation_days,1);
min_cons = zeros(1,num_days);
for i = 1:num_days
    % 1 p.m. to 5 p.m.
    idx_start = (i-1)*86400 + 1*3600;
    idx_stop = (i-1)*86400 + 5*3600;
    min_cons(i) = min(total_consumption(idx_start:idx_stop));
end
standby = median(min_cons) * length(total_consumption) / 3600;
inferred(length(experiments) + 3) = standby;
appliance_names_inferred = appliance_names;
appliance_names_inferred{length(experiments) + 3} = 'Standby';

%% Other
other_inferred = sum(total_consumption) / 3600 - sum(inferred);
inferred(length(experiments) + 4) = other_inferred;
other_actual = sum(total_consumption) / 3600 - sum(actual);
actual(length(experiments) + 3) = other_actual;
appliance_names_inferred{length(experiments)+4} = 'Other';
appliance_names{length(experiments)+3} = 'Other';

%% Change TV and Laptop for better plotting
tmp = appliance_names(7);
appliance_names(7) = appliance_names(8);
appliance_names(8) = tmp;
tmp = appliance_names_inferred(7);
appliance_names_inferred(7) = appliance_names_inferred(8);
appliance_names_inferred(8) = tmp;
tmp = actual(7);
actual(7) = actual(8);
actual(8) = tmp;
tmp = inferred(7);
inferred(7) = inferred(8);
inferred(8) = tmp;

%% Save as csv file
fid = fopen([plot_folder, 'actual.csv'], 'w');
fprintf(fid, 'label,consumption,share\n');
for i = 1:length(actual)
    fprintf(fid, '%s, %.2f, %.3f\n', appliance_names{i}, actual(i)/1000, actual(i) / sum(actual));
end
fclose(fid);
fid = fopen([plot_folder, 'inferred.csv'], 'w');
fprintf(fid, 'label,consumption,share\n');
for i = 1:length(inferred)
    fprintf(fid, '%s, %.2f, %.3f\n', appliance_names_inferred{i}, inferred(i)/1000, inferred(i) / sum(inferred));
end
fclose(fid);

%
% fig = figure;
% hold on;
% subplot(1,2,1);
% pie(actual);
% title('Ground truth');
% subplot(1,2,2);
% pie(inferred);
% labels = appliance_names_inferred;
% legend(labels);
% title('Inferred');
% fig = make_report_ready(fig, 'size', [width, height], 'fontsize', fontsize);
% % print('-depsc2', '-cmyk', '-r600', [plot_folder, filename, '.eps']); % if eps is needed
% saveas(fig, [plot_folder, plot_filename, '.png'], 'png');
% close(fig);

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
function f = dtlz4a(x)

alpha = 10;

%% g function
sum = 0;
for i = 3:8
sum = sum + (x(:,i)-0.5)^2;
end
g = sum;

%% cost functions
f(:,1) = (1+g)*cos(((x(:,1)^alpha)*pi)/2)*cos(((x(:,2)^alpha)*pi)/2);
f(:,2) = (1+g)*cos(((x(:,1)^alpha)*pi)/2)*sin(((x(:,2)^alpha)*pi)/2);
f(:,3) = (1+g)*sin(((x(:,1)^alpha)*pi)/2);
end

author: Eric-Bradford

repo: TS-EMO


text:
function DataMat = in_data_mat( DataFile )
% IN_DATA_MAT: Read EEG recordings stored in a free .MAT file.

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2009-2016

%% ===== GET OPTIONS =====
% Get options
OPTIONS = bst_get('ImportEegRawOptions');

%% ===== READ MATRIX FILE =====
% Read file
RawMat = load(DataFile);
% Get fields
fields = fieldnames(RawMat);
validFields = {};
% Loop to find the possible fields
for i = 1:length(fields)
    if ~isempty(RawMat.(fields{i})) && isnumeric(RawMat.(fields{i})) && (numel(RawMat.(fields{i})) > 31)
        validFields{end+1} = fields{i};
    end
end
if isempty(validFields)
    DataMat = [];
    bst_error(['No valid recordings field in: "' DataFile '"'], 'Import EEG data', 0);
    return
end
% Ask user which field to use, if there is more than one
if (length(validFields) > 1)
    res = java_dialog('question', 'Please select the field that contains your EEG recordings:', ...
                      'Import EEG data', [], validFields);
    % If user did not answer: exit
    if isempty(res)
        DataMat = [];
        return
    end
else
    res = validFields{1};
end
% Use selected field
FileData = RawMat.(res);

% Check matrix orientation
switch OPTIONS.MatrixOrientation
    case 'channelXtime'
        % OK
    case 'timeXchannel'
        % Transpose needed
        FileData = permute(FileData, [2 1 3]);
end

% Build time vector
Time = ((0:size(FileData,2)-1) ./ OPTIONS.SamplingRate - OPTIONS.BaselineDuration);
% ChannelFlag
ChannelFlag = ones(size(FileData,1), 1);
% Apply voltage units (in Brainstorm: recordings are stored in Volts)
switch (OPTIONS.VoltageUnits)
    case '\muV'
        FileData = FileData * 1e-6;
    case 'mV'
        FileData = FileData * 1e-3;
    case 'V'
        % Nothing to change
    case 'None'
        % Nothing to change
end
% If only one time frame: double it
if (size(FileData, 2) == 1)
    FileData = repmat(FileData, [1 2 1]);
end

% If loading 3D matrix: 3rd dimension is the epoch list
nbEpoch = size(FileData, 3);

% Initialize returned structure
DataMat = db_template('datamat');
DataMat.Comment     = 'EEG/MAT';
DataMat.ChannelFlag = ChannelFlag;
DataMat.Time        = Time;
DataMat.Device      = 'Unknown';
DataMat.nAvg        = OPTIONS.nAvg;
DataMat = repmat(DataMat, [nbEpoch, 1]);

% Process each epoch
BaseComment = DataMat(1).Comment;
for i = 1:nbEpoch
    DataMat(i).F = double(FileData(:,:,i));
    % Add indice number for multiple epochs
    if (nbEpoch > 1)
        DataMat(i).Comment = sprintf('%s #%d', BaseComment, i);
    end
end

author: brainstorm-tools


text:
%  Figure 10.57      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%
% fig10_57.m is a script for disk read/write head control design
% the design is based on 1/s^2 with robustness
% with respect to a resonance of damping zeta = z at wm.
% the design calls for a single lead selected to maximize wc
% the phase margin is set by alpha = a
% the resonance is 'gain stabilized' with GM
% the time is scaled to milliseconds
clf;
wm= 5*pi
z=.05;
GM = 4;
a=0.1;
numGo=1;
denGo= [1 0 0];
sysGo=tf(numGo,denGo);
numG= [1/(50*pi) 1];
denG=[1/(25*pi^2) 1/(50*pi) 1  0 0];
sysG=tf(numG,denG);
b=GM/(2*z);
% PM >50; select alpha = 0.1
T= (1/wm)*sqrt(b/a^1.5);
wc = 1/(sqrt(a)*T)
% K = sqrt(a)*wc^2
K= sqrt(a)*wm^2/b;
numD =K*sqrt(a)*[T 1];
denD = [a*T 1];
sysD=tf(numD,denD);
sysOLo=sysD*sysGo;
sysOL=sysD*sysG;
figure(1)
clf
% bode(sysOLo);
% hold on
w=logspace(-1,2);
[mag,ph]=bode(sysOL,w);
subplot(211);
% convert to db
magdb=20*log10(mag(:,:));
magdb1=[magdb; zeros(size(magdb))];
semilogx(w,magdb1);
grid;
title('Fig. 10.57 Disk drive control with gain-stabilized resonance');
xlabel('\omega (rad/sec)');
ylabel('Magnitude (db)');
subplot(212);
semilogx(w,[ph(:,:); -180*ones(size(ph(:,:)))]);
grid;
xlabel('\omega (rad/sec)');
ylabel('Phase (deg)');
pause;
%hold off
[Gm,Pm,Wcg,Wcp] = margin(mag,ph,w)

author: Sable

repo: mcbench-benchmarks


text:
function [ lam, mu ] = cycle_brent ( f, x0 )

%*****************************************************************************80
%
%% CYCLE_BRENT finds a cycle in an iterated mapping using Brent's method.
%
%  Discussion:
%
%    Suppose we a repeatedly apply a function f(), starting with the argument
%    x0, then f(x0), f(f(x0)) and so on.  Suppose that the range of f is finite.
%    Then eventually the iteration must reach a cycle.  Once the cycle is reached,
%    succeeding values stay within that cycle.
%
%    Starting at x0, there is a "nearest element" of the cycle, which is
%    reached after MU applications of f.
%
%    Once the cycle is entered, the cycle has a length LAM, which is the number
%    of steps required to first return to a given value.
%
%    This function uses Brent's method to determine the values of MU and LAM,
%    given F and X0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Richard Brent,
%    An improved Monte Carlo factorization algorithm,
%    BIT,
%    Volume 20, Number 2, 1980, pages 176-184.
%
%  Parameters:
%
%    Input, integer F ( integer I ), the name of the function
%    to be analyzed.
%
%    Input, integer X0, the starting point.
%
%    Output, integer LAM, the length of the cycle.
%
%    Output, integer MU, the index in the sequence starting
%    at X0, of the first appearance of an element of the cycle.
%
  power = 1;
  lam = 1;
  tortoise = x0;
  hare = f ( x0 );

  while ( tortoise ~= hare )
    if ( power == lam )
      tortoise = hare;
      power = power * 2;
      lam = 0;
    end
    hare = f ( hare );
    lam = lam + 1;
  end

  mu = 0;
  tortoise = x0;
  hare = x0;

  for i = 0 : lam - 1
    hare = f ( hare );
  end

  while ( tortoise ~= hare )
    tortoise = f ( tortoise );
    hare = f ( hare );
    mu = mu + 1;
  end

  return
end


text:
function dimensions=getPathwayDimensions(pathway)
% getPathwayDimensions
%   Retrieves the dimensions of metabolic network in a pathway structure.
%   Returns the position of the upper left corner, width and height.
%
%   pathway         pathway structure representing the pathway to be drawn
%
%   dimension       a 1x4 vector with x and y for the upper left corner,
%                   height and width
%
%   Usage: dimensions=getPathwayDimensions(pathway)

right=0;
left=inf;
top=inf;
bottom=0;

%Loops through the compartments to find the right and bottom border and the
%position of the upper left corner
for i=1:length(pathway.listOfCompartments)
    if pathway.listOfCompartments(1,i).x<left
        left=pathway.listOfCompartments(1,i).x;
    end
    if pathway.listOfCompartments(1,i).y<top
        top=pathway.listOfCompartments(1,i).y;
    end
    if (pathway.listOfCompartments(1,i).x+pathway.listOfCompartments(1,i).w)>right
        right=pathway.listOfCompartments(1,i).x+pathway.listOfCompartments(1,i).w;
    end
    if (pathway.listOfCompartments(1,i).y+pathway.listOfCompartments(1,i).h)>bottom
        bottom=pathway.listOfCompartments(1,i).y+pathway.listOfCompartments(1,i).h;
    end
end

%Loops through the species to find the object furthest to the right, left,
%top and bottom
for i=1:length(pathway.listOfSpecies)
    if pathway.listOfSpecies(1,i).x<left
        left=pathway.listOfSpecies(1,i).x;
    end
    if pathway.listOfSpecies(1,i).y<top
        top=pathway.listOfSpecies(1,i).y;
    end
    if (pathway.listOfSpecies(1,i).x+pathway.listOfSpecies(1,i).w)>right
        right=pathway.listOfSpecies(1,i).x+pathway.listOfSpecies(1,i).w;
    end
    if (pathway.listOfSpecies(1,i).y+pathway.listOfSpecies(1,i).h)>bottom
        bottom=pathway.listOfSpecies(1,i).y+pathway.listOfSpecies(1,i).h;
    end
end

dimensions=[left,top,right-left,bottom-top];
end

author: SysBioChalmers

repo: RAVEN


text:
mins=60;
debug.set_debug(1);

%% Sakurajima
%chanmatch = 'chan=~/[BESH]H[ENZ]/  || chan=~/BD[FL123]/';
%ds = datasource('antelope', ...
       %'/raid/data/sakurajima/Seismic_Infrasound/SAK_dbs/GT/dbSAK_GT');
%setup('Sakurajima', datenum(2015,5,18), datenum(2015,6,7), 'pf/setup_Sakurajima.pf', chanmatch)
%iceweb(ds, 'thissubnet', 'Sakurajima', 'snum', datenum(2015,5,18), 'enum', datenum(2015,6,7), 'delaymins', 0, 'matfile', 'pf/Sakurajima.mat', 'nummins', mins, 'runmode', 'archive');

%% Montserrat
% chanmatch = 'chan=~/[BESH]H[ENZ]/  || chan=~/BD[FL123]/';
%chanmatch = 'chan=~/[BESH]H[ENZ]/';
chanmatch = 'chan=~/[BESH]HZ/';
ds = datasource('antelope', ...
        '/raid/data/MONTSERRAT/antelope/db/db%04d%02d%02d',...
        'year','month','day');
setup('Montserrat', datenum(1995,7,1), datenum(2009,1,1), 'pf/setup_Montserrat.pf', chanmatch)
%montserrat_remove_analog_sites
iceweb(ds, 'thissubnet', 'Montserrat', 'snum', datenum(1998,2,19), 'enum', datenum(1998,2,19,23,59,59), 'delaymins', 0, 'matfile', 'pf/Montserrat.mat', 'nummins', mins, 'runmode', 'test');

%% Others yet to be converted to iceweb
% NOTE:
%   iceweb is just run in a single matlab session without waveform mat files and no rtexec
%   tremor_loadwaveformdata is for running concurrently with an rtexec system mixed with live data and creates waveform mat files)
%tremor_loadwaveformdata('thissubnet', 'Pavlof', 'snum', datenum(2013, 5, 17), 'enum', datenum(2013, 5, 18), 'delaymins', 0, 'matfile', 'pf/tremor_runtime.mat', 'nummins', mins, 'runmode', 'archive');
%tremor_loadwaveformdata('thissubnet', 'Pavlof', 'snum', datenum(2013, 5, 16), 'enum', datenum(2013, 5, 17), 'delaymins', 0, 'matfile', 'pf/tremor_runtime.mat', 'nummins', mins, 'runmode', 'archive');
%tremor_loadwaveformdata('thissubnet', 'Pavlof', 'snum', datenum(2013, 5, 15), 'enum', datenum(2013, 5, 16), 'delaymins', 0, 'matfile', 'pf/tremor_runtime.mat', 'nummins', mins, 'runmode', 'archive');
%tremor_loadwaveformdata('thissubnet', 'Pavlof', 'snum', datenum(2013, 5, 14), 'enum', datenum(2013, 5, 15), 'delaymins', 0, 'matfile', 'pf/tremor_runtime.mat', 'nummins', mins, 'runmode', 'archive');

%tremor_loadwaveformdata('thisubnet', 'Redoubt', 'snum', datenum(2009, 3, 23, 0, 0, 0), 'enum', datenum(2009, 4, 5), 'nummins', mins, 'runmode', 'archive');
%tremor_loadwaveformdata('thisubnet', 'Redoubt', 'snum', datenum(2009, 3, 19), 'enum', datenum(2009, 3, 23), 'nummins', mins, 'runmode', 'archive');
%tremor_loadwaveformdata('thisubnet', 'Redoubt', 'snum', datenum(2009, 1, 1, 0, 0, 0), 'enum', datenum(2009, 3, 19), 'nummins', mins, 'runmode', 'archive');
%tremor_loadwaveformdata('thisubnet', 'Redoubt', 'snum', datenum(2009, 4, 5), 'enum', datenum(2009, 4, 12), 'nummins', mins, 'runmode', 'archive');
%tremor_loadwaveformdata('thisubnet', 'Redoubt', 'snum', datenum(2009, 4, 12), 'enum', datenum(2009, 6, 1), 'nummins', mins, 'runmode', 'archive');
%tremor_loadwaveformdata('thisubnet', 'Redoubt', 'snum', datenum(2008, 11, 1), 'enum', datenum(2009, 1, 1), 'nummins', mins, 'runmode', 'archive');
%tremor_loadwaveformdata('thisubnet', 'Redoubt', 'snum', datenum(2008, 9, 1), 'enum', datenum(2008, 11, 1), 'nummins', mins, 'runmode', 'archive');
%tremor_loadwaveformdata('thisubnet', 'Okmok', 'snum', datenum(2008, 7, 15), 'enum', datenum(2008, 8, 20), 'nummins', mins, 'runmode', 'archive');
%tremor_loadwaveformdata('thisubnet', 'Pavlof', 'snum', datenum(2007, 8, 11), 'enum', datenum(2007, 9, 15), 'nummins', mins, 'runmode', 'archive');
%tremor_loadwaveformdata('thisubnet', 'Augustine', 'snum', datenum(2005, 12, 1), 'enum', datenum(2006, 3, 15), 'nummins', mins, 'runmode', 'archive');
%tremor_loadwaveformdata('thisubnet', 'Shishaldin', 'snum', datenum(1999, 2, 15), 'enum', datenum(1999, 5, 27), 'nummins', mins, 'runmode', 'archive');
%tremor_loadwaveformdata('thisubnet', 'Pavlof', 'snum', datenum(1996, 9, 10), 'enum', datenum(1997, 1, 1), 'nummins', mins, 'runmode', 'archive');

author: geoscience-community-codes

repo: GISMO


text:
function mu = trapezmf(z,a,b,c,d)
%TRAPEZMF Trapezoidal membership function.
%   MU = TRAPEZMF(Z,A,B,C,D) computes a fuzzy membership function with a
%   trapezoidal shape. Z is the input variable and can be a vector of
%   any length. A, B, C, and D are scalar parameters that define the
%   trapezoidal shape. The parameters must be ordered so that A <= B, B
%   <= C, and C <= D.
%
%       MU = 0                          Z < A
%       MU = (Z - A) ./ (B - A)         A <= Z < B
%       MU = 1                          B <= Z < C
%       MU = 1 - (Z - C) ./ (D - C)     C <= Z < D
%       MU = 0                          D <= Z
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

mu = zeros(size(z));

up_ramp_region = (a <= z) & (z < b);
top_region = (b <= z) & (z < c);
down_ramp_region = (c <= z) & (z < d);

mu(up_ramp_region) = 1 - (b - z(up_ramp_region))./(b - a);
mu(top_region) = 1;
mu(down_ramp_region) = 1 - (z(down_ramp_region) - c)./(d - c);

author: dipum

repo: dipum-toolbox

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


text:
clear;

% LPV model
lpv = {...
    @(p)1      @(p)p(1);
    @(p)p(2)^2 @(p)p(1)^2+p(2);
};
% number of states (size of the A matrix)
n = 4;
% parameter dependencies:
% dep(i,j,k) is 1 if Sp{i,j} depends on p(k)
dep = zeros([size(lpv) 2]);
dep(1,2,:) = [1 0];
dep(2,1,:) = [0 1];
dep(2,2,:) = [1 1];

% sampling intervals for each parameter
domain = [-5 5; -5 5];
% grid size: number of grid points for each parameter
gridsize = [11 11];

% hosvd_lpv
D1 = sampling_lpv(lpv, dep, domain, gridsize);
[S1 U1] = hosvd_lpv(D1, dep, gridsize);

% hosvd of a full array
f = @(p) [lpv{1,1}(p) lpv{1,2}(p); lpv{2,1}(p) lpv{2,2}(p)];
D2 = sampling_vec(f, domain, gridsize);
[S2 U2] = hosvd(D2, [1 1 0 0], 0.01);

% signs in U and S might be different
Sdiff = abs(S1(:)) - abs(S2(:));

% should be small
disp ' '
disp 'error (should be small):'
disp(norm(Sdiff))


text:
clear, close all
        % Determinarea numarului de noduri pe unitatea de lungime
nrnod=50;
        % Generarea unei retele pe domeniul considerat
[x,y]=meshgrid(0:1/(nrnod-1):2);
        % Determinarea numarului de noduri pe o linoe
N=length(x);
        % Initializarea matricelor V si Vp
V=zeros(N); Vp=V;
        % Fixarea erorii maxime impuse
erormax=0.0001;
        % Initializarea contorului de cicluri
nrp=1;
        % Ciclul repetitiv de calcul al potentialului
while nrp<3 | abs(max(max(V(2:N-1,2:N-1)-...
        Vp(2:N-1,2:N-1))))>erormax
        % - memorarea valorilor de la iteratia precedenta
    Vp=V;
        % - recalcularea potentialului in nodurile considerate
    V=V+del2(V);
	    % - impunerea conditiilor de frontiera
    V(N,:)=0; V(:,1)=0; V(:,N)=0; V(1,:)=1;
        % - incremenatarea contorului de cicluri
	nrp=nrp+1;
end;
        % Reprezentarea liniilor echipotentiale selectate
[c,h]=contour(x,y,V,[1,0.8,0.4,0.25,0.005]);
        % Etichetarea liniilor echipotentiale
clabel(c,h)
        % Setarea axelor patrate si personalizarea graficii
axis square
grid
xlabel('x')
ylabel('y')


text:
%   CLASS Remote_Resource_Manager
% =========================================================================
%
% DESCRIPTION
%
% EXAMPLE
%
% FOR A LIST OF CONSTANTs and METHODS use doc Remote_Resource_Manager

%--------------------------------------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:        Giulio Tagliaferro, Andrea Gatti ...
%  Contributors:      Giulio Tagliaferro, Andrea Gatti ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

classdef Remote_Resource_Manager < Ini_Manager
    properties (Constant, Access = private)
        DEFAULT_RESOURCE_FILE = 'remote_resource.ini';
        DEFAULT_RESOURCE_TXT = 'to be filled';
    end

    properties (Access = private)
        local_storage = '';
        credentials;
    end

    % =========================================================================
    %  INIT
    % =========================================================================
    methods (Access = public)
        % Guard the constructor against external invocation.  We only want
        % to allow a single instance of this class.  See description in
        % Singleton superclass.
        function this = Remote_Resource_Manager(file_name)
            if (nargin == 0)
                file_name = Remote_Resource_Manager.DEFAULT_RESOURCE_FILE;
                if exist(file_name, 'file') ~= 2
                    %Remote_Resource_Manager.writeDefault(); %deafult file not
                    %stable enough
                end
            end
            this = this@Ini_Manager(file_name);

            if ispc()
                home = [getenv('HOMEDRIVE') getenv('HOMEPATH')];
                this.local_storage = [home '\AppData\Local\goGPS'];
            else
                home = getenv('HOME');
                if ismac()
                    this.local_storage = [home '/Library/Application Support/goGPS'];
                else
                    this.local_storage = [home '/.goGPS'];
                end
            end
            if ~(exist(this.local_storage, 'dir'))
                mkdir(this.local_storage)
            end

            this.readFile();
            credentials_path = [Core.getInstallDir filesep 'credentials.txt'];
            if exist(credentials_path, 'file') == 0
                Core.getLogger.addError('The "credentials.txt" file is missing.\nIt will be created as empty from "credentials.example.txt" in goGPS folder');
                try
                    credentials_default_path = [Core.getInstallDir filesep 'credentials.example.txt'];
                    copyfile(credentials_default_path, credentials_path);
                    this.credentials = Ini_Manager(credentials_path);
                    this.credentials.readFile();
                catch
                    this.credentials = Ini_Manager();
                end
            else
                this.credentials = Ini_Manager(credentials_path);
                this.credentials.readFile();
            end
        end
    end

    % =========================================================================
    %  SINGLETON GETTERS
    % =========================================================================
    methods (Static)
        % Concrete implementation.  See Singleton superclass.
        function this = getInstance(resources_file, force_clean)
            % Get the persistent instance of the class
            persistent unique_instance_resource_manager__
            ini_is_present = false;
            switch nargin
                case 0
                    force_clean = false;
                case 1
                    if ischar(resources_file)
                        ini_is_present = true;
                        force_clean = false;
                    else
                        force_clean = resources_file;
                    end
                case 2
                    ini_is_present = true;
            end

            if force_clean
                log = Core.getLogger();
                clear unique_instance_resource_manager__;
                unique_instance_resource_manager__ = [];
            end

            if isempty(unique_instance_resource_manager__)
                if ini_is_present
                    this = Remote_Resource_Manager(resources_file);
                else
                    this = Remote_Resource_Manager();
                end
                unique_instance_resource_manager__ = this;
            else
                if ini_is_present
                    this = Remote_Resource_Manager(resources_file);
                    unique_instance_resource_manager__ = this;
                else
                    this = unique_instance_resource_manager__;
                end
            end
        end
    end

    methods
        function [ip, port, user, passwd] = getServerIp(this, name)
            % Return the ip of a server given the server name
            %
            % SYNTAX:
            %   [ip, port] = this.getServerIp(name)
            ip = [];
            port = [];
            ip_port = this.getData('SERVER', name);
            ip = ip_port{1};
            port = ip_port{2};
            this.credentials.readFile(); % read every time for updated credentials
            user = this.credentials.getData(name,'username');
            passwd = this.credentials.getData(name,'password');
        end

        function f_struct = getFileLoc(this, file_name)
            % Return the remote path of the file
            %
            % SYNTAX:
            %   f_struct = this.getFileLoc(file_name)
            f_struct.filename = this.getData(['f_' file_name],'filename');
            f_struct.const = this.getData(['f_' file_name],'sys');
            locations = this.getData(['f_' file_name],'location');
            if ~iscell(locations)
                locations = {locations};
            end
            f_struct.loc_number = length(locations);
            for i = 1 : f_struct.loc_number
                f_struct.(['loc' sprintf('%03d',i)]) = this.getData('LOCATION',locations{i});
            end

        end

        function [file_structure, latency] = getFileStr(this, center_name, resource_name)
            % Get the logical file structure for the desidered center and
            % resource the latecncy of the resource
            %
            % SYNTAX
            %   [file_structure, latency] = this.getFileStr(center_name, resource_name)
            %
            % OUTPUT:
            % file_strcuture
            %           the structure is a tree and can cointains fields
            %           named in 3 way 'fn' (where n is progressive number at
            %           the same level of the structure) , 'and' and 'or'.
            %           Or means that all sub field of the structure has to
            %           be found or means at least one.
            %           Leaves of the tree are cell containing the file code to
            %           be found in the remote resource ini file and a boolean to
            %           tell if the file has been found or not
            %
            %           example: fs.and.f1
            %                          .f2.or.f1
            %                                .f2
            %                          .f3
            %
            %                    f1 = {'cnes_erp' , 0}
            % latency
            %           [h1 h2] h1 -> hours before which we now the
            %                         resource is not there
            %                   h2 -> hours after which we are sure we will
            %                         found the resource
            %

            % Check if the requested resource is iono
            if iscell(center_name)
                center_name = center_name{1};
            end

            if numel(resource_name) > 5 && strcmp(resource_name(1:5), 'iono_')
                str = this.getData(['ic_' center_name], resource_name);
            elseif strcmp(resource_name, 'bias')
                str = this.getData(['bc_' center_name], resource_name);
            else
                str = this.getData(['oc_' center_name], resource_name);
            end
            if isempty(str)
                str = '';
                if ~strcmp(center_name, 'default')
                    str = ', using default';
                end
                if numel(resource_name) < 3 || ~strcmp(resource_name(1:3), 'vmf')
                    Core.getLogger.addWarning(sprintf('No resource %s for center %s%s',resource_name, center_name, str))
                end
                file_structure = [];
                latency = [];
            else
                file_structure = this.parseLogicTree(str);
                if numel(resource_name) > 5 && strcmp(resource_name(1:5), 'iono_')
                    latency = 0; % unkown latency (at the moment it is not present)
                elseif strcmp(resource_name, 'bias')
                    latency = this.getData(['bc_' center_name], [resource_name '_latency']);
                else
                    latency = this.getData(['oc_' center_name], [resource_name '_latency']);
                end
            end
        end
    end

    % Specific file query methods
    methods
        function [center, center_ss] = getCenterList(this, type)
            % Get the list of available centers and the supported constellations
            %
            % INPUT
            %   type    0: ORBIT (default)
            %           1: IONO
            %
            % SYNTAX
            %   [center, center_ss] = this.getCenterList(type)
            if nargin == 1 || isempty(type)
                type = 0;
            end
            if type == 1
                tmp = this.getData('IONO_CENTER', 'available');
            elseif type == 2
                tmp = this.getData('BIAS_CENTER', 'available');
            else
                tmp = this.getData('ORBIT_CENTER', 'available');
            end
            n_center = numel(tmp);
            center = cell(n_center, 1);
            center_ss = cell(n_center, 1);
            for c = 1 : n_center
                center{c} = regexp(tmp{c}, '(?<=@).*', 'match', 'once');
                if isempty(center{c})
                    center{c} = tmp{c};
                end
                center_ss{c} = regexp(tmp{c}, '.*(?=@)', 'match', 'once');
            end
        end

        function [center, center_ss] = getCenterListExtended(this, type)
            % Get the list of available centers with their description and the supported constellations
            %
            % INPUT
            %   type    0: ORBIT (default)
            %           1: IONO
            %
            % SYNTAX
            %   [center, center_ss] = this.getCenterList(type)
            if nargin == 1 || isempty(type)
                type = 0;
            end
            if type == 1
                tmp = this.getData('IONO_CENTER', 'available');
            elseif type == 2
                tmp = this.getData('BIAS_CENTER', 'available');
            else
                tmp = this.getData('ORBIT_CENTER', 'available');
            end
            n_center = numel(tmp);
            center = cell(n_center, 1);
            center_ss = cell(n_center, 1);
            for c = 1 : n_center
                if type == 1 || type == 2
                    center{c} = tmp{c};
                else
                    center{c} = regexp(tmp{c}, '(?<=@).*', 'match', 'once');
                end
                if isempty(center{c})
                    center{c} = tmp{c};
                else
                    if type == 1
                        center{c} = [upper(center{c}) ' - ' this.getData(['ic_' center{c}], 'description')];
                    elseif type == 2
                        center{c} = [upper(center{c}) ' - ' this.getData(['bc_' center{c}], 'description')];
                    else
                        center{c} = [upper(center{c}) ' - ' this.getData(['oc_' center{c}], 'description')];
                    end
                end
                center_ss{c} = regexp(tmp{c}, '.*(?=@)', 'match', 'once');
            end
        end

        function [descr] = centerToString(this, center)
            % Get the center description
            %
            % SYNTAX
            %   descr = this.centerToString(center)
            descr = sprintf('Current active center: "%s" - %s\n\nAvailable resources:\n', center, this.getData(['oc_' center], 'description'));

            % get resource_list
            key = this.getKeys(['oc_' center]);
            for k = 1 : numel(key)
                if isempty(regexp(key{k}, '(description)|(.*_latency)', 'once'))
                    if strcmp(key{k}, 'iono_center')
                        % get resource_list
                        iono_center = this.getData(['oc_' center], 'iono_center');
                        key_iono = this.getKeys(['ic_' iono_center]);
                        for i = 1 : numel(key_iono)
                            if isempty(regexp(key_iono{i}, '(description)|(.*_latency)', 'once'))
                                resource = this.resourceTreeToString(['ic_' iono_center], key_iono{i});
                                descr = sprintf('%s%s', descr, resource);
                            end
                        end
                    else
                        descr = sprintf('%s%s', descr, this.resourceTreeToString(['oc_' center], key{k}));
                    end
                end
            end

            if not(strcmp(center, 'default'))
                center = 'default';
                descr = sprintf('%s\n\n\nResources fallback: "%s" - %s\n\nAvailable resources:\n', descr, center, this.getData(['oc_' center], 'description'));

                % get resource_list
                key = this.getKeys(['oc_' center]);
                for k = 1 : numel(key)
                    if isempty(regexp(key{k}, '(description)|(.*_latency)', 'once'))
                        if strcmp(key{k}, 'iono_center')
                            % get resource_list
                            iono_center = this.getData(['oc_' center], 'iono_center');
                            key_iono = this.getKeys(['ic_' iono_center]);
                            for i = 1 : numel(key_iono)
                                if isempty(regexp(key_iono{i}, '(description)|(.*_latency)', 'once'))
                                    resource = this.resourceTreeToString(['ic_' iono_center], key_iono{i});
                                    descr = sprintf('%s%s', descr, resource);
                                end
                            end
                        else
                            descr = sprintf('%s%s', descr, this.resourceTreeToString(['oc_' center], key{k}));
                        end
                    end
                end
            end
        end

        function [flag_frub] = getOrbitType(this, center)
            % Get the orbit type availability
            %
            % SYNTAX
            %   [flag_frub] = this.getOrbitType(center)
            if iscell(center)
                center = center{1};
            end
            flag_frub(1) = ~isempty(this.getData(['oc_' center], 'final'));
            flag_frub(2) = ~isempty(this.getData(['oc_' center], 'rapid'));
            flag_frub(3) = ~isempty(this.getData(['oc_' center], 'ultra'));
            flag_frub(4) = ~isempty(this.getData(['oc_' center], 'broadcast'));
            flag_frub(5) = ~isempty(this.getData(['oc_' center], 'real-time'));
        end

        function [flag_fp1p2b] = getIonoType(this, iono_center)
            % Get the iono type availability
            %
            % SYNTAX

            %   [flag_frub] = this.getIonoType(center)
            state = Core.getCurrentSettings();
            if nargin == 1
                iono_center = state.getRemoteIonoCenter();
            end
            if isempty(iono_center)
                iono_center = 'default';
            end

            flag_fp1p2b = false(5, 1);

            flag_fp1p2b(1) = ~isempty(this.getData(['ic_' iono_center], 'iono_final')) && ~strcmp(this.getData(['ic_' iono_center], 'iono_final'), 'empty');
            flag_fp1p2b(2) = ~isempty(this.getData(['ic_' iono_center], 'iono_rapid')) && ~strcmp(this.getData(['ic_' iono_center], 'iono_rapid'), 'empty');
            flag_fp1p2b(3) = ~isempty(this.getData(['ic_' iono_center], 'iono_predicted1')) && ~strcmp(this.getData(['ic_' iono_center], 'iono_predicted1'), 'empty');
            flag_fp1p2b(4) = ~isempty(this.getData(['ic_' iono_center], 'iono_predicted2')) && ~strcmp(this.getData(['ic_' iono_center], 'iono_predicted2'), 'empty');
            flag_fp1p2b(5) = ~isempty(this.getData(['ic_' iono_center], 'iono_broadcast')) && ~strcmp(this.getData(['ic_' iono_center], 'iono_broadcast'), 'empty');

            if ~any(flag_fp1p2b) && ~strcmp(iono_center, 'none')
                % Switch to default center
                flag_fp1p2b(1) = ~isempty(this.getData(['ic_default'], 'iono_final'));
                flag_fp1p2b(2) = ~isempty(this.getData(['ic_default'], 'iono_rapid'));
                flag_fp1p2b(3) = ~isempty(this.getData(['ic_default'], 'iono_predicted1'));
                flag_fp1p2b(4) = ~isempty(this.getData(['ic_default'], 'iono_predicted2'));
                flag_fp1p2b(5) = ~isempty(this.getData(['ic_default'], 'iono_broadcast'));
            end
        end

        function [flag_frub] = getVMFResType(this)
            % Get the vmf res type availability
            %
            % SYNTAX
            %   [flag_frub] = this.getVMFResType(center)

            state = Core.getCurrentSettings();

            if state.mapping_function == 1
                flag_frub = [false false false];
            elseif state.mapping_function == 2
                flag_frub = [false true false];
            elseif state.mapping_function == 3
                flag_frub = [false false false];
            elseif state.mapping_function == 4
                flag_frub = [true false true];
            end
        end

        function [flag_frub] = getVMFSourceType(this)
            % Get the vmf availability type availability
            %
            % SYNTAX
            %   [flag_frub] = this.getVMFSourceType(center)
            state = Core.getCurrentSettings();

            if state.mapping_function == 1
                flag_frub = [false false false];
            elseif state.mapping_function == 2
                flag_frub = [true false true];
            elseif state.mapping_function == 3
                flag_frub = [false false false];
            elseif state.mapping_function == 4
                flag_frub = [true true true];
            elseif state.mapping_function == 5
                flag_frub = [true true true];
            end
        end

        function [tree_str] = resourceTreeToString(this, center, resource_type)
            % Parse a resource tree
            %
            % SYNTAX
            %   [tree_str] = this.resourceTreeToString(center, resource_type)
            tmp = this.getData(center, resource_type);
            tree = this.parseLogicTree(tmp);
            if iscell(tree)
                clear tmp
                tmp.f1 = tree;
                tree = tmp;
            end
            tree_str = sprintf(' - %s\n', resource_type);
            if ~isempty(tree)
                tree_str = sprintf('%s%s', tree_str, this.treeToString(tree, 4));
            end
        end

        function remote_path = resourceToString(this, resource)
            if strcmp(resource, 'null')
                remote_path = 'none';
            else
                try
                    file_name = this.getData(['f_', resource], 'filename');
                    location = this.getData(['f_', resource], 'location');
                    if iscell(location)
                        location = location{1};
                    end
                    if isempty(location)
                        location = '';
                        server = cell(2,1);
                    else
                        location = this.getData('LOCATION', location);
                        % remove the server
                        server = regexp(location, '(?<=\?\{)[a-z\_A-Z]*(?=\})', 'match');
                        if isempty(server)
                            server = cell(2,1);
                        else
                            location = regexp(location, '(?<=(\?\{[a-z_A-Z]*\})).*', 'match');
                        end
                        server = this.getData('SERVER', server{1});
                    end
                    % The protocol could be improved woth more values
                    if iscell(server) && numel(server) == 2 && ~isempty(server{2})
                        switch server{2}
                            case '21'
                                protocol = 'ftp://';
                            otherwise
                                protocol = 'http://';
                        end
                    else
                        protocol = '';
                    end
                    remote_path = [protocol server{1} ':' server{2} location{1} file_name];
                catch
                    remote_path = sprintf('corrupted resource "%s" -> check resource file', resource);
                end
            end
        end

        function tree_str = treeToString(this, tree, indentation_lev)
            % Parse a resource tree
            %
            % SYNTAX
            %   [tree_str] = this.treeToString(tree, indentation_lev)
            if nargin == 1
                indentation_lev = 0;
            end
            str_padding = char(32*ones(1, indentation_lev + 1));
            tree_str = '';
            tree_names = fieldnames(tree);
            for b = 1 : numel(tree_names)
                if this.isNode(tree_names(b))
                    % is a node
                    tree_str = sprintf('%s%s|- %s\n%s', tree_str, str_padding, upper(tree_names{b}), this.treeToString(tree.(tree_names{b}), indentation_lev + 4));
                else
                    % is a leaf
                    if iscell(tree.(tree_names{b}))
                        % tree_str = sprintf('%s%s|- %s\n', tree_str, str_padding, tree.(tree_names{b}){1});
                        tree_str = sprintf('%s%s|- %s\n', tree_str, str_padding, this.resourceToString(tree.(tree_names{b}){1}));
                    elseif isstruct(tree.(tree_names{b}))
                        tree_str = sprintf('%s%s', tree_str, this.treeToString(tree.(tree_names{b}), indentation_lev));
                    else
                        % NON VALID
                    end
                end
            end
        end
    end

    methods ( Access = public)
        function file_structure = parseLogicTree(this, str)
            % Description parse the logic structure found in
            % remote_resource.ini to get the file structure descripbed in
            % this.getFileStr
            %
            % SYNTAX:
            %   file_structure = this.parseLogicTree(this)
            [status, list] = this.findElements(str);
            file_structure = [];
            if status == 0
                file_structure  = {strtrim(str), false};
                return
            else
                if status == 1
                    cond = 'or';
                else
                    cond = 'and';
                end
                for i = 1 : length(list)
                    file_structure.(cond).(['f' num2str(i)]) = this.parseLogicTree(list{i});
                end
            end

        end

        function [status, list] = findElements(this, str)
            % Return the element of the string separated by | and &
            %
            % SYNTAX:
            %       [status, list] = this.findElements(str)
            %
            % OUTPUT:
            % status: -1 and 0 nothing 1 or
            % list: list of string parts
            [matches] = regexp(str, '\&|\|', 'match');
            if isempty(matches)
                status = 0;
                list = {};
                return
            else
                open = 0; %number if open pharentesis
                status = '';
                index = [];
                for i = 1:length(str)
                    if open == 0
                        if str(i) == '&' | str(i) == '|'
                            if isempty(status)
                                status = str(i);
                                index = [index; i];
                            else
                                if status ~= str(i)
                                    Core.getLogger.addWarning('| and & can not exist at the same level, check parenthesis')
                                    status = 0;
                                    return
                                else
                                    index = [index; i];
                                end
                            end
                        end
                    end
                    if str(i) == '('
                        open = open + 1;
                    end
                    if str(i) == ')'
                        open = open - 1;
                    end
                end
                if status == '|'
                    status = 1;
                elseif status == '&'
                    status = -1;
                end
                list = {};
                for i = 1 : length(index)
                    if i == 1
                        list{end + 1} = this.removeTrailingPar(str(1:index(i)-1));
                    else
                        list{end + 1} = this.removeTrailingPar(str(index(i-1)+1: index(i)-1));
                    end
                end
                list{end + 1} = this.removeTrailingPar(str(index(end)+1 : end));

            end
        end

        function str = removeTrailingPar(this, str)
            % Remove trailing parenthesis
            %
            % SYNTAX:
            %   str = this.removeTrailingPar(str)
            for i =1 :length(str)
                if str(i)~=' '
                    if str(i)=='('
                        str(1:i) = [];
                    end
                    break
                end
            end
            for i =length(str) : -1:1
                if str(i)~=' '
                    if str(i)==')'
                        str(i:end) = [];
                    end
                    break
                end
            end
        end
    end

    methods (Static)
        function writeDefault(this)
            % Write the deafut remote resource ini file if it is not found
            %
            % SYNTAX:
            %       Remote_Resource_Manager.writeDefault()
            fid = fopen(Remote_Resource_Manager.DEFAULT_RESOURCE_FILE,'w+');
            fprintf(fid, Remote_Resource_Manager.DEFAULT_RESOURCE_TXT);
            fclose(fid);
        end

        function is_node = isNode(tree_el)
            % Local usage return if a node of node_tree is an "and" or an "or"
            if numel(tree_el) == 1 && iscell(tree_el) && ~isempty(regexp(tree_el{1}, '(and)*|(or)*'))
                is_node = true;
            else
                is_node = false;
            end
        end

    end
end

author: goGPS-Project


text:
function flg = checkDictUse
% Since MATLAB 7.3(2006A) and 7.4(2006B) doesnot support dicom dictionary
% this function checks the matlab version and returns the flag to allow the
% use of dictionary.
%
% DK
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

% set flag to 1 so as to always use Dictionary
flg = 1;

mathVer = version;

if str2num(mathVer(1:3)) > 7.3
    flg = 0;
end

return

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
mcc -m computeFeaturesCompiled ...
  -R -singleCompThread ...
  -a filehandling/ ...
  -a misc/ ...
  -a toolbox/ ...
  -a optflow_deqing/ ...
  -a deepmatching/deepmex_cluster.mexa64
if exist('readme.txt','file'),
  delete('readme.txt');
end
if exist('mccExcludedFiles.log','file'),
  delete('mccExcludedFiles.log');
end
if exist('requiredMCRProducts.txt','file'),
  delete('requiredMCRProducts.txt');
end

author: kristinbranson

repo: JAABA


text:
function moon_lander()
close all;

t0 = 0;
tfinal = 800;
initial_state = [0 160000 0 0 -100 0 15000]';

[t, x] = runge_kutta(@moonlander, initial_state, [t0 tfinal], 0.05);

figure, plot(x(1,:), x(2,:), '-b*')
title('Landing trajectory')

figure,
plot(t, x(3,:), 'k')
title('Attitude')

figure,
plot(t, x(5,:), 'r')
title('Landing speed')

figure,
plot(t, x(7,:), 'g')
title('Mass')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Helper function to solve a second order differential equation.
% Called from function exerciseB()
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function xd = moonlander(t, x)
%target speed and x position
vtarget = -100;
xtarget = 0;

%constants
mass_min = 8000;
g = 1.6;
J = 100000;
ce=3.0*1000;
thlmax = 1000*0.5;
thtmax = 1000*44;
% compute control actions!
Kt = 300;
Kl = -0.000001;
vy = x(5);
posx = x(1);
posy = x(2);

if posy<3000
ut = Kt*(vtarget-vy);
ul = Kl*(xtarget-posx);
ul = sat(ul, -1, 1);
ut = sat(ut, 0, 1);
else
    ut=0;
    ul=0;
end

if x(7)<=mass_min
    Fl = 0;
    Ft = 0;
else
    Fl=thlmax*ul;
    Ft=thtmax*ut;
end

xd1 = x(4);
xd2 = x(5);
xd3 = x(6);
xd4 = (1/x(7))*(Fl*cos(x(3)) - Ft*sin(x(3)));
xd5 = (1/x(7))*(Fl*sin(x(3)) + Ft*cos(x(3))) - g;
xd6 = 4*Fl/J;

if x(7)<=mass_min
    % not changing mass
    display('ran out of fuel')
    % the derivative of mass is zero
    xd7 = 0 ;
else
    xd7 = -(abs(Fl/ce)+abs(Ft/ce));
end
xd = [xd1 xd2 xd3 xd4 xd5 xd6 xd7]';
return

function u= sat(u, min_val, max_val)
if u> max_val
    u = max_val;
end
if u < min_val
    u = min_val;
end

repo: ARTE


text: a=input('Introduceti valoarea lui a:')


text:
function [document, scores] = report_precision_recall(context, experiment, trackers, sequences, varargin)
% report_overlap Generate a report using tracking precision recall methodology
%
% Performs tracking precision-recall analysis and generates a report based on the results.
%
% Input:
% - context (structure): Report context structure.
% - experiment (struct): An experiment structure.
% - trackers (cell): An array of tracker structures.
% - sequences (cell): An array of sequence structures.
% - varargin[UseTags] (boolean): Analyze according to tags (otherwise according to sequences).
% - varargin[HideLegend] (boolean): Hide legend in plots.
%
% Output:
% - document (structure): Resulting document structure.
% - scores (struct): A scores structure.
%

usetags = get_global_variable('report_tags', true);
hidelegend = get_global_variable('report_lagend_hide', false);
resolution = 100;

for i = 1:2:length(varargin)
    switch lower(varargin{i})
        case 'usetags'
            usetags = varargin{i+1};
        case 'resolutuion'
            resolution = varargin{i+1};
        case 'hidelegend'
            hidelegend = varargin{i+1};
        otherwise
            error(['Unknown switch ', varargin{i}, '!']) ;
    end
end

if ~strcmp(experiment.type, 'unsupervised')
   error('Tracking precision-recall analysis only suitable for unsupervised experiments!');
end

document = document_create(context, 'tpr', 'title', 'Tracking precision recall');

trackers_hash = md5hash(strjoin((cellfun(@(x) x.identifier, trackers, 'UniformOutput', false)), '-'), 'Char', 'hex');
parameters_hash = md5hash(sprintf('%d%d', usetags, resolution));

tags = {};

if isempty(experiment.tags)
    usetags = false;
end;

if usetags && isfield(experiment, 'tags')
    tags = union(experiment.tags, {'all'});
    sequences_hash = md5hash(strjoin(tags, '-'), 'Char', 'hex');
else
    sequences_hash = md5hash(strjoin((cellfun(@(x) x.name, sequences, 'UniformOutput', false)), '-'), 'Char', 'hex');
end;

cache_identifier = sprintf('tpr_%s_%s_%s_%s', experiment.name, trackers_hash, sequences_hash, parameters_hash);

result = document_cache(context, cache_identifier, @analyze_precision_recall, experiment, trackers, ...
    sequences, 'Tags', tags, 'Resolution', resolution);

if usetags
    % When using tags we have inserted a separate one for this
    mask = strcmp('tag_all', result.selectors);

    average_curve = result.curves(:, mask);
    average_measures = result.measures(:, mask);

    % Now remove the 'all' tag from results
    tag_curve = result.curves(:, ~mask);
    tag_measures = result.measures(:, ~mask);

    selector_tags = cat(2, result.selectors(~mask), result.selectors(mask));

else

    average_curve = cell(numel(trackers), 1);
    average_measures = zeros(numel(trackers), 3);

    for t = 1:numel(trackers)
        average_curve{t} = mean(cat(3, result.curves{t, :}), 3);
        f = 2 * (average_curve{t}(:, 1) .* average_curve{t}(:, 2)) ./ (average_curve{t}(:, 1) + average_curve{t}(:, 2));
        [average_measures(t, 1), idx] = max(f);
        average_measures(t, 2) = average_curve{t}(idx, 1);
        average_measures(t, 3) = average_curve{t}(idx, 2);
    end;

    tag_curve = result.curves;
    tag_measures = result.measures;

    selector_tags = result.selectors;

end

scores.name = 'TPR';
scores.values = average_measures;
scores.ids = {'f', 'tp', 'tr'};
scores.names = {'F', 'TP', 'TR'};
scores.order = {'descending', 'descending', 'descending'};

tracker_labels = cellfun(@(x) iff(isfield(x.metadata, 'verified') && x.metadata.verified, [x.label, '*'], x.label), trackers, 'UniformOutput', 0);

print_text('Writing tracking precision-recall table ...');

document.section('Experiment %s', experiment.name);

pr_plot(document, sprintf('%s_average', experiment.name), ...
    sprintf('Experiment %s (average)', experiment.name), ...
    trackers, average_curve, hidelegend);

table_data = highlight_best_rows(num2cell(cat(2, tag_measures(:, : , 1), average_measures(:, 1))), repmat({'descending'}, 1, size(tag_measures, 2) + 1));

document.table(table_data, 'columnLabels', selector_tags, 'rowLabels', tracker_labels, 'title', 'Tracking precision-recall overview');

document.subsection('Detailed plots');

for t = 1:size(tag_curve, 2)

    plot_title = sprintf('Tracking precision-recall plot for tag %s in experiment %s', ...
        selector_tags{t}, experiment.name);
    plot_id = sprintf('overlap_%s_%s', experiment.name, selector_tags{t});

    pr_plot(document, plot_id, plot_title, trackers, tag_curve(:, t), ~hidelegend);

end;

document.write();

end

function pr_plot(document, identifier, title, trackers, curves, hidelegend)

    handle = plot_blank('Visible', false, 'Title', 'Overlap', 'Width', 6, 'Height', 6); hold on;

    phandles = zeros(numel(trackers), 1);

    for t = 1:numel(curves)
        phandles(t) = plot(curves{t}(:, 2), curves{t}(:, 1), 'Color', trackers{t}.style.color);
    end;

    labels = cellfun(@(x) x.label, trackers, 'UniformOutput', false);

    if ~hidelegend
        legend(phandles, labels);
    end;

    xlabel('Tracking recall');
    ylabel('Tracking precision');
    xlim([0, 1]);
    ylim([0, 1]);
    hold off;
    document.figure(handle, identifier, title);

    close(handle);
end

author: votchallenge

repo: toolkit-legacy


text:
function task_division_test ( )

%*****************************************************************************80
%
%% TASK_DIVISION_TEST tests the TASK_DIVISION library.
%
%  Discussion:
%
%    This program simply demonstrates how one might automate the
%    assignment of T tasks to P processors, assuming that the assignment
%    is to be beforehand.
%
%    In that case, we just want to make sure that we assign each task
%    to a processor, that we assign about the same number of tasks
%    to each processor, and that we assign each processor a contiguous
%    range of tasks, say tasks I_LO to I_HI.
%
%    The routine that is called simulates this process.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 October 2011
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TASK_DIVISION_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the TASK_DIVISION library.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Demonstrate how to automate the division of\n' );
  fprintf ( 1, '  T tasks among a range of P processors\n' );
  fprintf ( 1, '  indexed from PROC_FIRST to PROC_LAST.\n' );

  task_number = 23;
  proc_first = 0;
  proc_last = 3;
  task_division ( task_number, proc_first, proc_last );

  task_number = 17;
  proc_first = 1;
  proc_last = 6;
  task_division ( task_number, proc_first, proc_last );

  task_number = 17;
  proc_first = 4;
  proc_last = 6;
  task_division ( task_number, proc_first, proc_last );

  task_number = 5;
  proc_first = -2;
  proc_last = 6;
  task_division ( task_number, proc_first, proc_last );

  task_number = 5;
  proc_first = 0;
  proc_last = 4;
  task_division ( task_number, proc_first, proc_last );

  task_number = 5;
  proc_first = 0;
  proc_last = 0;
  task_division ( task_number, proc_first, proc_last );

  task_number = 1000;
  proc_first = 1;
  proc_last = 17;
  task_division ( task_number, proc_first, proc_last );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TASK_DIVISION_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function obj = JZMO_wbmethod(obj)

global V_ropt W_ropt Nk Nrf Ns n Vn H;
t1 = clock;
[V_RF, V_U] =  MO_AltMinWB(V_ropt);
[W_RF, W_B] =  MO_AltMinWB(W_ropt);
%W_B = W_B/20;
v = zeros(Nk,1);
V_B = zeros(Nrf,Ns,Nk);
for i = 1:Nk
    v(i) = trace(V_U(:,:,i)'*V_RF'*V_RF*V_U(:,:,i));
    V_B(:,:,i)= V_U(:,:,i) /sqrt(v(i));
end
for k = 1:Nk
    newH(:,:,k) = W_B(:,:,k)'*W_RF'*H(:,:,k)*V_RF*V_B(:,:,k);
    W_X(:,:,k) = inv(newH(:,:,k)*newH(:,:,k)' +  Vn * v(k) *W_B(:,:,k)'*(W_RF)'*W_RF*W_B(:,:,k))*newH(:,:,k);
    newnewH(:,:,k) = W_X(:,:,k)'*newH(:,:,k);
    W_B(:,:,k) = W_B(:,:,k) * W_X(:,:,k) ;
end
t2 = clock;
runtime  = etime(t2,t1);
obj.V_B = V_B;
obj.W_B = W_B;
obj.V_RF = V_RF;
obj.W_RF = W_RF;
obj.runtime = obj.runtime + runtime;
obj = get_wbmetric(obj);
% if (obj.ber(n)>30)
%     save H
%     pause
% end
a = 1;

author: Zzhaoxingyu

repo: hybrid-beamforming-for-three-scenes


text:
function headmodel = ft_headmodel_singlesphere(mesh, varargin)

% FT_HEADMODEL_SINGLESPHERE creates a volume conduction model of the
% head by fitting a spherical model to a set of points that describe
% the head surface.
%
% For MEG this implements Cuffin BN, Cohen D.  "Magnetic fields of
% a dipole in special volume conductor shapes" IEEE Trans Biomed Eng.
% 1977 Jul;24(4):372-81.
%
% Use as
%   headmodel = ft_headmodel_singlesphere(mesh, ...)
%
% Optional arguments should be specified in key-value pairs and can include
%   conductivity     = number, conductivity of the sphere
%
% See also FT_PREPARE_VOL_SENS, FT_COMPUTE_LEADFIELD

% FIXME document both EEG and MEG case

% Copyright (C) 2012-2013, Donders Centre for Cognitive Neuroimaging, Nijmegen, NL
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% get the optional arguments
conductivity = ft_getopt(varargin, 'conductivity', 1);

if any(strcmp(varargin(1:2:end), 'unit')) || any(strcmp(varargin(1:2:end), 'units'))
  % the geometrical units should be specified in the input mesh
  ft_error('the ''unit'' option is not supported any more');
end

if isnumeric(mesh) && size(mesh,2)==3
  % assume that it is a Nx3 array with vertices
  % convert it to a structure, this is needed to determine the units further down
  mesh = struct('pos', mesh);
elseif isstruct(mesh) && isfield(mesh,'bnd')
  % take the triangulated surfaces from the input structure
  mesh = mesh.bnd;
end

% replace pnt with pos
mesh = fixpos(mesh);

if ~isstruct(mesh) || numel(mesh)>1 || ~isfield(mesh, 'pos')
  ft_error('the input mesh should be a set of points or a single triangulated surface')
end

if numel(conductivity)~=1
  ft_error('the conductivity should be a single number')
end

if numel(mesh)~=1
  ft_error('fitting a single sphere requires a single mesh')
end

% start with an empty volume conductor
headmodel = [];

% ensure that the mesh has units, estimate them if needed
mesh = ft_determine_units(mesh);

% copy the geometrical units into the volume conductor
headmodel.unit = mesh.unit;

% fit a single sphere to all headshape points
[single_o, single_r] = fitsphere(mesh.pos);

headmodel.r    = single_r;
headmodel.o    = single_o;
headmodel.cond = conductivity;
headmodel.type = 'singlesphere';

fprintf('single sphere: radius = %.1f, conductivity = %f\n', headmodel.r, headmodel.cond);

author: spm


text:
function failed_bug2418

% WALLTIME 00:10:00
% MEM 2gb
% DEPENDENCY ft_senstype ft_senslabel

%% test the consistency between labels and type

type = {
  'ant128'
  'biosemi64'
  'biosemi128'
  'biosemi256'
  'bti148'
  'bti148_planar'
  'bti248'
  'bti248_planar'
%  'btiref'         % returns a subset
  'ctf151'
  'ctf151_planar'
  'ctf275'
  'ctf275_planar'
%  'ctfheadloc'     % returns a subset
%  'ctfref'         % returns a subset
  'eeg1005'
  'eeg1010'
  'eeg1020'
%  'ext1020'        % returns lower and upper case
  'egi32'
  'egi64'
  'egi128'
  'egi256'
  'neuromag122'
%  'neuromag122alt' % does not apply any more, see bugzilla
  'neuromag306'
%  'neuromag306alt' % does not apply any more, see bugzilla
  'itab28'
  'itab153'
  'itab153_planar'
  'yokogawa9'
  'yokogawa64'
  'yokogawa64_planar'
  'yokogawa160'
  'yokogawa160_planar'
  'yokogawa440'
  };

for i=1:length(type)
  lab = ft_senslabel(type{i});
  lab = lab(:);

  disp(type{i});
  assert(strcmp(type{i}, ft_senstype(lab)));
end

%% test the provided problematic data

cd(dccnpath('/home/common/matlab/fieldtrip/data/test'));
load bug2418.mat

assert(isequal(ft_senstype(testdata_short), 'neuromag306'));
assert(isequal(ft_senstype(testdata_short.label), 'neuromag306'));
assert(isequal(ft_senstype(testdata_short.grad), 'neuromag306'));
assert(isequal(ft_senstype(testdata_short.grad.label), 'neuromag306'));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ frames, indexes ] = vec2frames( vec, Nw, Ns, direction, window, padding )
% VEC2FRAMES Splits signal into overlapped frames using indexing.
%
%   B=vec2frames(A,M,N) creates a matrix B whose columns consist of
%   segments of length M, taken at every N samples along input vector A.
%
%   [B,R]=vec2frames(A,M,N,D,W,P) creates a matrix B whose columns
%   or rows, as specified by D, consist of segments of length M, taken
%   at every N samples along the input vector A and windowed using the
%   analysis window specified by W. The division of A into frames is
%   achieved using indexes returned in R as follows: B=A(R);
%
%   Summary
%
%           A is an input vector
%
%           M is a frame length (in samples)
%
%           N is a frame shift (in samples)
%
%           D specifies if the frames in B are rows or columns,
%             i.e., D = 'rows' or 'cols', respectively
%
%           W is an optional analysis window function to be applied to
%             each frame, given as a function handle, e.g., W = @hanning
%             or as a vector of window samples, e.g., W = hanning( M )
%
%           P specifies if last frame should be padded to full length,
%             or simply discarded, i.e., P = true or false, respectively
%
%           B is the output matrix of frames
%
%           R is a matrix of indexes used for framing, such that division
%             of A into frames is achieved as follows: B=A(R);
%
%   Examples
%
%           % divide the input vector into seven-sample-long frames with a shift
%           % of three samples and return frames as columns of the output matrix
%           % (note that the last sample of the input vector is discarded)
%           vec2frames( [1:20], 7, 3 )
%
%           % divide the input vector into seven-sample-long frames with a shift
%           % of three samples and return frames as rows of the output matrix
%           % (note that the last sample of the input vector is discarded)
%           vec2frames( [1:20], 7, 3, 'rows' )
%
%           % divide the input vector into seven-sample-long frames with a shift
%           % of three samples, pad the last frame with zeros so that no samples
%           % are discarded and return frames as rows of the output matrix
%           vec2frames( [1:20], 7, 3, 'rows', [], true )
%
%           % divide the input vector into seven-sample-long frames with a shift
%           % of three samples, pad the last frame with white Gaussian noise
%           % of variance (1E-5)^2 so that no samples are discarded and
%           % return frames as rows of the output matrix
%           vec2frames( [1:20], 7, 3, 'rows', false, { 'noise', 1E-5 } )
%
%           % divide the input vector into seven-sample-long frames with a shift
%           % of three samples, pad the last frame with zeros so that no samples
%           % are discarded, apply the Hanning analysis window to each frame and
%           % return frames as columns of the output matrix
%           vec2frames( [1:20], 7, 3, 'cols', @hanning, 0 )
%
%   See also FRAMES2VEC, DEMO

%   Author: Kamil Wojcicki, UTD, July 2011

    % usage information
    usage = 'usage: [ frames, indexes ] = vec2frames( vector, frame_length, frame_shift, direction, window, padding );';

    % default settings
    switch( nargin )
    case { 0, 1, 2 }, error( usage );
    case 3, padding=false; window=false; direction='cols';
    case 4, padding=false; window=false;
    case 5, padding=false;
    end

    % input validation
    if( isempty(vec) || isempty(Nw) || isempty(Ns) ), error( usage ); end;
    if( min(size(vec))~=1 ), error( usage ); end;
    if( Nw==0 || Ns==0 ), error( usage ); end;

    vec = vec(:);                       % ensure column vector

    L = length( vec );                  % length of the input vector
    M = floor((L-Nw)/Ns+1);             % number of frames

    % perform signal padding to enable exact division of signal samples into frames
    % (note that if padding is disabled, some samples may be discarded)
    if( ~isempty(padding) )

        % figure out if the input vector can be divided into frames exactly
        E = (L-((M-1)*Ns+Nw));

        % see if padding is actually needed
        if( E>0 )

            % how much padding will be needed to complete the last frame?
            P = Nw-E;

            % pad with zeros
            if( islogical(padding) && padding )
                vec = [ vec; zeros(P,1) ];

            % pad with a specific numeric constant
            elseif( isnumeric(padding) && length(padding)==1 )
                vec = [ vec; padding*ones(P,1) ];

            % pad with a low variance white Gaussian noise
            elseif( isstr(padding) && strcmp(padding,'noise') )
                vec = [ vec; 1E-6*randn(P,1) ];

            % pad with a specific variance white Gaussian noise
            elseif( iscell(padding) && strcmp(padding{1},'noise') )
                if( length(padding)>1 ), scale = padding{2};
                else, scale = 1E-6; end;
                vec = [ vec; scale*randn(P,1) ];

            % if not padding required, decrement frame count
            % (not a very elegant solution)
            else
                M = M-1;

            end

            % increment the frame count
            M = M+1;
        end
    end

    % compute index matrix
    switch( direction )

    case 'rows'                                                 % for frames as rows
        indf = Ns*[ 0:(M-1) ].';                                % indexes for frames
        inds = [ 1:Nw ];                                        % indexes for samples
        indexes = indf(:,ones(1,Nw)) + inds(ones(M,1),:);       % combined framing indexes

    case 'cols'                                                 % for frames as columns
        indf = Ns*[ 0:(M-1) ];                                  % indexes for frames
        inds = [ 1:Nw ].';                                      % indexes for samples
        indexes = indf(ones(Nw,1),:) + inds(:,ones(1,M));       % combined framing indexes

    otherwise
        error( sprintf('Direction: %s not supported!\n', direction) );

    end

    % divide the input signal into frames using indexing
    frames = vec( indexes );

    % return if custom analysis windowing was not requested
    if( isempty(window) || ( islogical(window) && ~window ) ), return; end;

    % if analysis window function handle was specified, generate window samples
    if( isa(window,'function_handle') )
        window = window( Nw );
    end

    % make sure analysis window is numeric and of correct length, otherwise return
    if( isnumeric(window) && length(window)==Nw )

        % apply analysis windowing beyond the implicit rectangular window function
        switch( direction )
        case 'rows', frames = frames * diag( window );
        case 'cols', frames = diag( window ) * frames;
        end

    end

% EOF

author: a-nagrani

repo: VGGVox

save_path: github-repos/MATLAB/a-nagrani-VGGVox


text:
% END - last index to memmapdata array
%
% Author: Arnaud Delorme, SCCN, INC, UCSD, Nov. 2008

% Copyright (C) 2008 Arnaud Delorme, SCCN, INC, UCSD
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function s = end(a, k, n);

    s = size(a, k);

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function memory_test ( n_log_min, n_log_max )

%*****************************************************************************80
%
%% MEMORY_TEST carries out some tests of memory limits in MATLAB.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 June 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N_LOG_MIN, the logarithm of the minimum N to try.
%
%    Input, integer N_LOG_MAX, the logarithm of the maximum N to try.
%
  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'MEMORY_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Try to see how big vectors and matrices can be.\n' );

  if ( nargin < 1 )
    n_log_min = 0;
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Using default value of N_LOG_MIN = %d\n', n_log_min );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, '  User input value of N_LOG_MIN = %d\n', n_log_min );
  end

  if ( nargin < 2 )
    n_log_max = 25;
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Using default value of N_LOG_MAX = %d\n', n_log_max );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, '  User input value of N_LOG_MAX = %d\n', n_log_max );
  end
%
%  I4VEC test.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4VEC Memory Test\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, ...
    'Log2(N)            N   Ave     CPU        Real\n' );
  fprintf ( 1, '\n' );

  for n_log = n_log_min : n_log_max

    i4vec_memory_test ( n_log );

  end
%
%  R8VEC test.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8VEC Memory Test\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, ...
    'Log2(N)            N   Ave     CPU        Real\n' );
  fprintf ( 1, '\n' );

  for n_log = n_log_min : n_log_max

    r8vec_memory_test ( n_log );

  end
%
%  I4MAT test.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4MAT Memory Test\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, ...
    'Log2(N)            N            N1            N2   Ave     CPU        Real\n' );
  fprintf ( 1, '\n' );

  for n_log = n_log_min : n_log_max

    i4mat_memory_test ( n_log );

  end
%
%  R8MAT test.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8MAT Memory Test\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, ...
    'Log2(N)            N            N1            N2   Ave     CPU        Real\n' );
  fprintf ( 1, '\n' );

  for n_log = n_log_min : n_log_max

    r8mat_memory_test ( n_log );

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'MEMORY_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function i4mat_memory_test ( n_log )

%*****************************************************************************80
%
%% I4MAT_MEMORY_TEST declares and uses an I4MAT of total size N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 June 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N_LOG, the logarithm base 2 of N.
%
  n = 2^n_log;

  n1_log = floor ( n_log / 2 );
  n1 = 2^n1_log;
  n2_log = n_log - n1_log;
  n2 = 2^n2_log;

  fprintf ( 1, '  %4d  %12d', n_log, n );

  tic;
  cpu_time1 = cputime;

  i4mat = zeros ( n1, n2 );

  i4mat = round ( 2.0 * rand ( n1, n2 ) );

  average = ( sum ( sum ( i4mat(1:n1,1:n2) ) ) ) / n1 / n2;

  cpu_time2 = cputime;
  real_diff = toc;

  cpu_diff = cpu_time2 - cpu_time1;

  fprintf ( 1, '  %12d  %12d  %4.2f  %10e  %10e\n', ...
    n1, n2, average, cpu_diff, real_diff );

  return
end
function i4vec_memory_test ( n_log )

%*****************************************************************************80
%
%% I4VEC_MEMORY_TEST declares and uses an I4VEC of size N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 June 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N_LOG, the logarithm base 2 of N.
%
  n = 2^n_log;

  fprintf ( 1, '  %4d  %12d', n_log, n );

  tic;
  cpu_time1 = cputime;

  i4vec = zeros ( n, 1 );

  i4vec = round ( 2.0 * rand ( n, 1 ) );

  average = ( sum ( i4vec(1:n) ) ) / n;

  cpu_time2 = cputime;
  real_diff = toc;

  cpu_diff = cpu_time2 - cpu_time1;

  fprintf ( 1, '  %4.2f  %10e  %10e\n', average, cpu_diff, real_diff );

  return
end
function r8mat_memory_test ( n_log )

%*****************************************************************************80
%
%% R8MAT_MEMORY_TEST declares and uses an R8MAT of total size N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 June 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N_LOG, the logarithm base 2 of N.
%
  n = 2^n_log;

  n1_log = floor ( n_log / 2 );
  n1 = 2^n1_log;
  n2_log = n_log - n1_log;
  n2 = 2^n2_log;

  fprintf ( 1, '  %4d  %12d', n_log, n );

  tic;
  cpu_time1 = cputime;

  r8mat = zeros ( n1, n2 );

  r8mat = 2.0 * rand ( n1, n2 );

  average = ( sum ( sum ( r8mat(1:n1,1:n2) ) ) ) / n1 / n2;

  cpu_time2 = cputime;
  real_diff = toc;

  cpu_diff = cpu_time2 - cpu_time1;

  fprintf ( 1, '  %12d  %12d  %4.2f  %10e  %10e\n', ...
    n1, n2, average, cpu_diff, real_diff );

  return
end
function r8vec_memory_test ( n_log )

%*****************************************************************************80
%
%% R8VEC_MEMORY_TEST declares and uses an R8VEC of size N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 June 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N_LOG, the logarithm base 2 of N.
%
  n = 2^n_log;

  fprintf ( 1, '  %4d  %12d', n_log, n );

  tic;
  cpu_time1 = cputime;

  r8vec = zeros ( n, 1 );

  r8vec = 2.0 * rand ( n, 1 );

  average = ( sum ( r8vec(1:n) ) ) / n;

  cpu_time2 = cputime;
  real_diff = toc;

  cpu_diff = cpu_time2 - cpu_time1;

  fprintf ( 1, '  %4.2f  %10e  %10e\n', average, cpu_diff, real_diff );

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% test example
% April 15th, 2010, By Reza Farrahi Moghaddam, Synchromedia Lab, ETS, Montreal, Canada

%
u0_filename = 'P03.tif';

% method name. Please try one of these methods: 'min_avg'; % 'avg'; 'normal';
% method = 'min_avg' is recommended for document image processing with less dependency on color.
method_name = 'min_avg'; % 'min_avg'; % 'avg'; 'normal';

u0_color = imread(u0_filename);

% ref image, using the standard method.
u0_gray_ref = double(rgb2gray(u0_color));

% main function
u0_gray = universal_color_to_gray_converter(u0_color , 'method', method_name);

% display
figure, imshow(u0_gray_ref / max(max(u0_gray_ref)))
figure, imshow(double(u0_gray) / max(max(double(u0_gray))))

% save the output
imwrite(u0_gray, [u0_filename,'-',method_name,'.tif']);