text:
%% This script simulates the LODCO-based epsilon-Greedy algorithm.
%  author: Hailiang Zhao
clc, clear
opt = optimset('Display', 'none');

%% basic parameter settings (had better not change those paras)
k = 1e-28;                % effective switched capacitance (a constant decided by the chip architecture)
tau = 0.002;              % the length of time slot (in second)
phi = 0.002;              % the cost of task dropping (in second)
omega = 1e6;              % the bandwidth of MEC server (in Hz)
sigma = 1e-13;            % the noise power of the receiver (in W)
p_tx_max = 1;             % the maximum transmit power of mobile device (in W)
f_max = 1.5e9;            % the maximum CPU-cycle frequency of mobile device (in Hz)
E_max = 0.002;            % the maximum amout of battery output energy (in J)
L = 1000;                 % the input size of the computation task (in bit)
X = 737.5;                % the number of CPU cycles needed on processing one bit of task
W = L * X;                % the number of CPU cycles needed on processing one task
E_H_max = 48e-6;          % the upper bound of the energy arrive at the mobile device (in J)
p_H = E_H_max / (2*tau);  % the average Energy Harvesting (EH) power (in W)
g0 = power(10, -4);       % the path-loss constant
d0 = 1;                   % the relative distance between each mobile device and each MEC server

%% parameter control
N = 10;                   % the number of mobile devices
M = 8;                    % the number of MEC servers
T = 2000;                 % the number of time slot (a.k.a. the size of the time horizon)
tau_d = 0.002;            % execution deadline (in second)
d = 50;                   % the distance between the mobile device and the MEC server (in meter)
E_min = 0.02e-3;          % the minimum amout of battery output energy (in J)
V = 1e-5;                 % the weight of penalty (the control parameter introduced by Lyapunov Optimization)
rho = 0.6;                % the probability that the computation task is requested
max_connects = 4;         % the maximum number of processible mobile devices for each MEC server ($ \frac{f_s^{max} \tau}{L X} $)
min_distance = 10;        % the minimum distance from mobile device to MEC server
max_distance = 50;        % the maximum distance from mobile device to MEC server
eps = 0.8;                % if rand() <= eps, for the best i-j pair, we always choose MEC server execution mode (except no MEC server to choose)

% the lower bound of perturbation parameter
E_max_hat = min(max(k * W * (f_max)^2, p_tx_max * tau), E_max);
theta = E_max_hat + V * phi / E_min;

%% allocate storage for valuable results
B = zeros(T, N);                   % the battery energy level (in J)
B_hat = zeros(T, N);               % the virtual battery energy level ($B_hat = B - theta$)
e = zeros(T, N);                   % the amout of the harvested and stored energy (in J)
chosen_mode = zeros(T, N);         % {1: local, 2: remote, 3: drop, 4: no task request}
chosen_server = zeros(T, N);       % record the index of chosen server for each mobile device if its choice is MEC server execution
f = zeros(T, N);                   % the CPU-cycle frequency of local execution (in Hz)
p = zeros(T, N);                   % the transmit power of computation offloading (in W)

mobile_exe_cost = zeros(T, N);     % the mobile execution cost (delay) (in second)
server_exe_cost = zeros(T, N);     % the MEC server execution cost (delay) (in second)
final_chosen_cost = zeros(T, N);   % the final execution delay under currently chosen modes (in second)

mobile_exe_E = zeros(T, N);        % the energy consumption for mobile execution (in J)
server_exe_E = zeros(T, N);        % the energy consumption for MEC server execution (in J)
final_chosen_E = zeros(T, N);      % the energy consumption of the final chosen modes (in J)

%% simulation begin
t = 1;
while t <= T
    disp(['===> Time slot #', num2str(t), ' <==='])

    %% allocate storage for mode-chosen
    device_server_pairs = [];      % each column represents i, j, J_s^{\star}(i, j), respectively
    remained_connects = max_connects * ones(M, 1);    % the available connections of MEC servers
    J_m = zeros(N, 1); J_s = zeros(N, M);             % the matrices for J_m and J_s values
    p_mat = zeros(N, M);                              % the matrix for transmit power
    server_cost_mat = zeros(N, M);                    % the matrix for MEC server execution cost
    server_E_mat = zeros(N, M);                       % the matrix for energy consumption of MEC server

    %% initialization
    % generate the virtual battery energy level
    B_hat(t, :) = B(t, :) - theta;
    % generate the channel power gain (from each mobile device to each MEC sever)
    distances = unifrnd(min_distance, max_distance, N, M);
    gamma = exprnd(1, N, M);
    h_mat = g0 * gamma .* power(d0 ./ distances, 4);

    %% step 1: for each mobile device, choose the initial mode
    for i = 1: N
        disp(['Mobile device #', num2str(i)])

        %% step 1.1: get the optimal energy harvesting no matter whether task is requested
        E_H_t = unifrnd(0, E_H_max);
        if B_hat(t, i) <= 0
            e(t, i) = E_H_t;
        end

        %% step 1.2: get the (initial) optimal computation offloading strategy (I_m(i), I_s(i, :), I_d(i), f(t, i), p(t, i))
        % generate the task request
        zeta = binornd(1, rho);
        if zeta == 0
            % chosen mode has to be 4
            disp('no task request generated!')
            chosen_mode(t, i) = 4;
        else
            % chosen_mode is chosen from {1, 2, 3}
            disp('task request generated!')

            %% step 1.2.1: solve the optimization problem $\mathcal{P}_{ME}$ (f(t, i) > 0)
            % calculate f_L and f_U
            f_L = max(sqrt(E_min / (k * W)), W / tau_d);
            f_U = min(sqrt(E_max / (k * W)), f_max);
            if f_L <= f_U
                % the sub-problem is feasible
                disp('mobile execution ($\mathcal{P}_{ME}$) is feasible!')

                if B_hat(t, i) < 0
                    f_0 = (V / (-2 * B_hat(t, i) * k))^(1/3);
                else
                    % complex number may exist, which may lead to error
                    f_0 = -(V / (2 * B_hat(t, i) * k))^(1/3);
                end

                if (f_0 > f_U && B_hat(t, i) < 0) || (B_hat(t, i) >= 0)
                    f(t, i) = f_U;
                elseif f_0 >= f_L && f_0 <= f_U && B_hat(t, i) < 0
                    f(t, i) = f_0;
                elseif f_0 < f_L && B_hat(t, i) < 0
                    f(t, i) = f_L;
                end
                % check whether f(t, i) is zero
                if f(t, i) == 0
                    disp('Something wrong! f is 0!')
                end

                % calculate the delay of mobile execution
                mobile_exe_cost(t, i) = W / f(t, i);
                % calculate the energy consumption of mobile execution
                mobile_exe_E(t, i) = k * W * (f(t, i)^2);
                % calculate the value of optimization goal
                J_m(i) = -B_hat(t, i) * k * W * (f(t, i)^2 + V * W / f(t, i));
            else
                % the sub-problem is not fasible because (i) the limited
                % computation capacity or (ii) time cosumed out of deadline
                % or (iii) the energy consumed out of battery energy level
                % If it is not feasible, it just means that we cannot choose
                % 'I_m(i)=1'. It dosen't mean that the task has to be dropped.
                disp('mobile execution ($\mathcal{P}_{ME}$) is not feasible!')
                f(t, i) = 0;
                mobile_exe_cost(t, i) = 0;
                mobile_exe_E(t, i) = 0;
                % 'I_m(i)=1' can never be chosen if mobile execution goal is inf
                J_m(i) = inf;
            end

            %% step 1.2.2: solve the optimization problem $\mathcal{P}_{SE}$ (p(t, i) > 0)
            % calculate J_s(i, j) from mobile device i to each MEC server j
            for j = 1: M
                disp(['MEC server #', num2str(j)])
                h = h_mat(i, j);

                E_tmp = sigma * L * log(2) / (omega * h);
                p_L_taud = (power(2, L / (omega * tau_d)) - 1) * sigma / h;
                % calculate p_L
                if E_tmp >= E_min
                    p_L = p_L_taud;
                else
                    % calculate p_E_min (use inline function and fsolve)
                    y = @(x) x * L - omega * log2(1 + h*x/sigma) * E_min;
                    % accroding to the function figure, p_L_taud is a positive
                    % number around 0.2
                    p_E_min = fsolve(y, 0.2, opt);
                    p_L = max(p_L_taud, p_E_min);
                end
                % calculate p_U
                if E_tmp >= E_max
                    p_U = 0;
                else
                    % caculate p_E_max (use inline function and fsolve)
                    y = @(x) x * L - omega * log2(1 + h*x/sigma) * E_max;
                    % accroding to the function figure, p_E_max is a large positive
                    % number around 20
                    p_E_max = fsolve(y, 100, opt);
                    p_U = min(p_tx_max, p_E_max);
                end

                if p_L <= p_U
                    % the sub-problem is feasible
                    disp('MEC server execution ($\mathcal{P}_{SE}$) is feasible!')
                    % calculate p_0
                    virtual_battery = B_hat(t, i);
                    y = @(x) virtual_battery * log2(1 + h*x/sigma) + ...
                        h * (V - virtual_battery*x) / log(2) / (sigma + h*x);
                    p_0 = fsolve(y, 0.5, opt);

                    if (p_U < p_0 && B_hat(t, i) < 0) || B_hat(t, i) >= 0
                        p_mat(i, j) = p_U;
                    elseif p_0 < p_L && B_hat(t) < 0
                        p_mat(i, j) = p_L;
                    elseif p_0 >= p_L && p_0 <= p_U && B_hat(t, i) < 0
                        p_mat(i, j) = p_0;
                    end
                    % check whether p_mat(i, j) is zero
                    if p_mat(i, j) == 0
                        disp('Something wrong! p is 0!')
                    end

                    % calculate the delay of MEC server execution
                    server_cost_mat(i, j) = L / (omega * log2(1 + h*p_mat(i, j)/sigma));
                    % calculate the energy consumption of MEC server execution
                    server_E_mat(i, j) = p_mat(i, j) * server_cost_mat(i, j);
                    % calculate the value of optimization goal
                    J_s(i, j) = (-B_hat(t, i) * p_mat(i, j) + V) * server_cost_mat(i, j);

                    % (we can not set server_exe_cost(t, i) and server_exe_E(t, i) for now)
                else
                    % the sub-problem is not feasible because (i) the limited transmit
                    % power or (ii) time cosumed out of deadline or (iii) the energy
                    % consumed out of battery energy level
                    % If it is not feasible, it just means that we cannot choose
                    % 'I_s(i,j)=1'. It dosen't mean that the task has to be dropped.
                    disp('MEC server execution ($\mathcal{P}_{SE}$) is not feasible!')
                    p_mat(i, j) = 0;
                    server_cost_mat(i, j) = 0;
                    server_E_mat(i, j) = 0;
                    % 'I_s(i,j)=1' can never be chosen if MEC server execution goal is inf
                    J_s(i, j) = inf;

                    % (we can not set server_exe_cost(t, i) and server_exe_E(t, i) for now)
                % Similarly, we do not check whether the energy cunsumed is larger than
                % battery energy level because the problem $\mathcal{J}_{CO}$ does
                % not have constraint (8).
                end
            end

            %% step 1.2.3: choose the (initial) optimal execution mode
            J_d = V * phi;
            disp(['J_m(i):', num2str(J_m(i))])
            disp(['J_s(i,:):', num2str(J_s(i, :))])
            [~, mode] = min([J_m(i), J_s(i, :), J_d]);
            if mode == 1
                % mobile execution
                chosen_mode(t, i) = 1;
                final_chosen_cost(t, i) = mobile_exe_cost(t, i);
                final_chosen_E(t, i) = mobile_exe_E(t, i);
            elseif mode == (M+2)
                % drop
                chosen_mode(t, i) = 3;
                final_chosen_cost(t, i) = phi;
                final_chosen_E(t, i) = 0;
            else
                % MEC servre execution
                chosen_mode(t, i) = 2;
                % add i, the chosen j, and their J_s(i, j) value into the pairs
                device_server_pairs = [device_server_pairs; [i, mode-1, J_s(i, mode-1)]];
            end
        end
    end

    %% step 2: allocate connection right for each mobile device who chooses MEC server execution
    while ~isempty(device_server_pairs)
        if rand() <= eps
            % find the {i, j, J_s(i, j)} with the minimum J_s(i, j) in the device_server_pairs, choose MEC server execution mode,
            % except no MEC server to choose
            [min_Js, idx] = min(device_server_pairs(:, 3));
            i = device_server_pairs(idx, 1); j = device_server_pairs(idx, 2);    % find the i and its j
            if remained_connects(j) >= 1
                % the chosen i can be allocated to the j-th MEC server
                % set its final modes as 2 and record the final cost and energy consumption for it
                chosen_mode(t, i) = 2;          % this is not necessary
                p(t, i) = p_mat(i, j);
                server_exe_cost(t, i) = server_cost_mat(i, j);
                server_exe_E(t, i) = server_E_mat(i, j);
                chosen_server(t, i) = j;
                final_chosen_cost(t, i) = server_exe_cost(t, i);
                final_chosen_E(t, i) = server_exe_E(t, i);

                % update remained_connects for j
                remained_connects(j) = remained_connects(j) - 1;
                % finally, remove it from global pairs, it is done
                device_server_pairs(device_server_pairs(:, 1) == i, :) = [];
            else
                % the chosen i can not be allocated to the j-th MEC server
                % set the mobile device's J_s(i, j) as inf
                % (which means no matter what happens, the MEC server j can never be chosen)
                J_s(i, j) = inf;
                % remove it from global pairs
                device_server_pairs(device_server_pairs(:, 1) == i, :) = [];
                % if there still exists j' where J_s(i, j') is not inf, we choose the best j',
                % otherwise we can only choose from mobile execution and drop
                if min(J_s(i, :)) ~= inf
                    [second_min_Js, second_min_j] = min(J_s(i, :));
                    % add the new pair into global pairs
                    device_server_pairs = [device_server_pairs; [i, second_min_j, second_min_Js]];
                else
                    % choose from mobile execution and drop for i
                    [~, mode] = min([J_m(i), inf, J_d]);
                    chosen_mode(t, i) = mode;
                    if mode == 1
                        final_chosen_cost(t, i) = mobile_exe_cost(t, i);
                        final_chosen_E(t, i) = mobile_exe_E(t, i);
                    else
                        final_chosen_cost(t, i) = phi;
                        final_chosen_E(t, i) = 0;
                    end
                end
            end
        else
            % do as LODCO-based Greedy algorithm do
            for j = 1: M
                % find every pair who chooses the MEC server j
                device_j_pairs = device_server_pairs(device_server_pairs(:, 2) == j, :);
                % find those devices is
                is = device_j_pairs(:, 1);
                if isempty(is)
                    disp(['For current MEC server #', num2str(j), ', no mobile device choose it!'])
                    % go to handle next MEC server
                    continue;
                end
                if remained_connects(j) >= length(is)
                    % every mobile device who chooses j can be connected,
                    % set their final modes as 2 and record the final cost and energy consumption for them
                    chosen_mode(t, is) = 2;     % this is not necessary
                    p(t, is) = transp(p_mat(is, j));
                    server_exe_cost(t, is) = transp(server_cost_mat(is, j));
                    server_exe_E(t, is) = transp(server_E_mat(is, j));
                    chosen_server(t, is) = repmat(j, 1, length(is));
                    final_chosen_cost(t, is) = server_exe_cost(t, is);
                    final_chosen_E(t, is) = server_exe_E(t, is);

                    % update remained_connects for j
                    remained_connects(j) = remained_connects(j) - length(is);
                    % finally, remove them from global pairs, they are done
                    device_server_pairs(device_server_pairs(:, 2) == j, :) = [];
                else
                    if remained_connects(j) == 0
                        % no mobile device can be connected, remove all mobile devices who chooses j from global pairs
                        device_server_pairs(device_server_pairs(:, 2) == j, :) = [];
                        % set those mobile devices' J_s(i, j) as inf
                        % (which means no matter what happens, the MEC server j can never be chosen)
                        J_s(is, j) = inf;
                        % choose mode again for those i and insert new potential pairs into global pairs
                        for idx = 1: numel(is)
                            i = is(idx);
                            [~, mode] = min([J_m(i), J_s(i, :), J_d]);
                            if mode == 1
                                chosen_mode(t, i) = 1;
                                final_chosen_cost(t, i) = mobile_exe_cost(t, i);
                                final_chosen_E(t, i) = mobile_exe_E(t, i);
                            elseif mode == (M+2)
                                chosen_mode(t, i) = 3;
                                final_chosen_cost(t, i) = phi;
                                final_chosen_E(t, i) = 0;
                            else
                                chosen_mode(t, i) = 2;
                                device_server_pairs = [device_server_pairs; [i, mode-1, J_s(i, mode-1)]];
                            end
                        end
                    else
                        % some mobile devices can be connected, set their final modes as 2 and remove them from global pairs
                        % besides, record the final cost and energy consumption for them
                        % for the left i', remove them from global pairs and set J_s(i', j) as inf,
                        % choose mode again for those i' and insert new potential pairs into global pairs

                        % sort the J_s(is, j) and return those lucky idxs
                        [~, idxs] = sort(device_j_pairs(:, 3));
                        for idx = 1: remained_connects(j)
                            i = idxs(idx);
                            chosen_mode(t, i) = 2;
                            p(t, i) = p_mat(i, j);
                            server_exe_cost(t, i) = server_cost_mat(i, j);
                            server_exe_E(t, i) = server_E_mat(i, j);
                            chosen_server(t, i) = j;
                            final_chosen_cost(t, i) = server_exe_cost(t, i);
                            final_chosen_E(t, i) = server_exe_E(t, i);

                            % remove i from global pairs
                            device_server_pairs(device_server_pairs(:, 1) == i, :) = [];
                        end

                        % update remained_connects for j
                        remained_connects(j) = 0;

                        % for those unlucky mobile devices i', set J_s(i', j) as inf (i' are in currently device_server_pairs(: , j) now)
                        residual_is = device_server_pairs(device_server_pairs(:, 2) == j, 1);
                        J_s(residual_is, j) = inf;
                        for idx = 1: numel(residual_is)
                            residual_i = residual_is(idx);
                            [~, mode] = min([J_m(residual_i), J_s(residual_i, :), J_d]);
                            if mode == 1
                                chosen_mode(t, residual_i) = 1;
                                final_chosen_cost(t, residual_i) = mobile_exe_cost(t, residual_i);
                                final_chosen_E(t, residual_i) = mobile_exe_E(t, residual_i);
                            elseif mode == (M+2)
                                chosen_mode(t, residual_i) = 3;
                                final_chosen_cost(t, residual_i) = phi;
                                final_chosen_E(t, residual_i) = 0;
                            else
                                chosen_mode(t, residual_i) = 2;
                                device_server_pairs = [device_server_pairs; [residual_i, mode-1, J_s(residual_i, mode-1)]];
                            end
                        end
                    end
                end
            end
        end
    end

    %% step 3: update the battery energy level and go to the next time slot
    B(t + 1, :) = B(t, :) - final_chosen_E(t, :) + e(t, :);
    t = t + 1;

end

%% step 4: evaluate the simulation results
% 1. the battery energy level vs. time slot
figure
plot(1:T, B(1:T, :))
hold on
plot(1:T, repmat(theta + E_H_max, [T, 1]), '-')
title('Envolution of battery energy level')
xlabel('time slot')
ylabel('battery energy level $B_t$ of each mobile device', 'Interpreter','latex')

% 2. the average execution cost vs. time slot
accumulated = 0;
average_cost = zeros(T, N);
figure
for i = 1: N
    % draw for each mobile device
    request_num = 0;
    for t = 1: T
        accumulated = accumulated + final_chosen_cost(t, i);
        if chosen_mode(t, i) ~= 4
            % there exists task request
            request_num = request_num + 1;
        end
        average_cost(t, i) = accumulated / request_num;
    end
    plot(1:T, average_cost(:, i));
    hold on
end
title('Envolution of average execution cost')
xlabel('time slot')
ylabel('average execution cost $\frac{1}{T} \sum_{t=0}^{T-1} cost^t$ of each mobile device', 'Interpreter','latex')

% 3. the average ratio of each chosen mode of the ith mobile device vs. time slot
average_ratio = zeros(T, 3);
mobile_exe = 0; server_exe = 0; drop = 0;
request_num = 0;
i = 1;          % we simply choose the first mobile device
for t = 1: T
    if final_chosen_cost(t, i) == 0
        continue
    else
        request_num = request_num + 1;
        if chosen_mode(t, i) == 1
            mobile_exe = mobile_exe + 1;
        elseif chosen_mode(t, i) == 2
            server_exe = server_exe + 1;
        else
            drop = drop + 1;
        end
    end
    average_ratio(t, :) = [mobile_exe, server_exe, drop] / request_num;
end
figure
plot(1:T, average_ratio(:, 1));
hold on
plot(1:T, average_ratio(:, 2));
hold on
plot(1:T, average_ratio(:, 3));
legend('mobile execution', 'MEC server execution', 'drop')
title('Envolution of average ratio of chosen modes')
xlabel('time slot')
ylabel('average  ratio of chosen modes $\frac{1}{T} \sum_{t=0}^{T-1} \{I_m^t, I_s^t, I_d^t\}$ of the i-th mobile device', 'Interpreter','latex')

author: hliangzhao

repo: Edge-Computing-Codes


text:
function [model, ll] = gpsimMapUpdateF(model, options)

% GPSIMMAPUPDATEF Update posterior mean of f.
% FORMAT
% DESC updates the stored version of the posterior mean for f in
% the GPSIMMAP model.
% ARG model : the model for which f is to be updated.
% ARG options : options vector.
% RETURN model : the model with f updated.
% RETURN ll : the log likelihood after update.
%
% SEEALSO : gpsimMapFunctionalLogLikeGradient,
% gpsimMapFunctionalLogLikeHessian, gpsimMapCreate
%
% COPYRIGHT : Neil D. Lawrence, 2006
%
% MODIFIED : Pei Gao, 2008

% SHEFFIELDML

display = options(1);
tolf = options(2);
tol = options(3);
iters = options(14);
if iters == 0
    iters = 100;
end
f = gpsimMapFunctionalExtractParam(model);

if isfield(model, 'priorProtein') && ~isempty(model.priorProtein)
  model.consLambda = 5;
  maxLambda = 8;
  nCons = length(model.priorProtein);
else
  model.consLambda = 0;
  maxLambda = 1;
end

llold = gpsimMapFunctionalLogLikelihood(model);

for j = 1:maxLambda
  %fprintf('lambda = %d; display = %d\n', lambda, display);
  %if model.isConcave
  flag = 0;
  for i = 1:iters
    ll = llold - 1;
    gf = gpsimMapFunctionalLogLikeGradients(model);
    if model.consLambda
      consMat = zeros(size(model.invCovf));
      for k=1:nCons
        ftimeIndex = find((model.priorProteinTimes(k)-model.mapt)==0);
        consMat(ftimeIndex, ftimeIndex) = model.consLambda;
      end
      hf = inv(model.invCovf+consMat);
    else
      hf = model.covf;
    end

    new_drn = (hf*gf')';
    factor = 1;
    count = 0;
    fold = f;
    while ll < llold
        f = fold + factor*new_drn;
        model = gpsimMapFunctionalExpandParam(model, f);
        ll = gpsimMapFunctionalLogLikelihood(model);
        lldiff = ll - llold;
        count = count + 1;
        if count > 0 & display
          fprintf('gpsimMapUpdateF, lldiff: %2.4f, factor %2.4f\n', lldiff, ...
                  factor);
        end
        factor = factor/2;
        if lldiff < tol & max(abs(fold - f))<tolf
            flag = 1;
            break
        end
    end
    if display
      fprintf('Iteration %d, log likelihood %2.4f\n', i, ll);
    end
    llold = ll;
    if flag == 1
      break
    end
  end
  model.consLambda = model.consLambda*4;
end

if isempty(find(model.f))
  warning('TF is zero!');
end
%else
%  error('Not yet implemented non-concave update')
%end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%% QUAD_JOB_LOCAL divides the quadrature problem into tasks.
%
%  Discussion:
%
%    Each task is given a portion of the interval [0,1] over which
%    to integrate.
%
%    Each task returns the estimated integral over its subinterval.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 March 2010
%
%  Author:
%
%    John Burkardt
%
  clear

  fprintf ( 1, '\n' );
  fprintf ( 1, 'QUAD_JOB_LOCAL\n' );
  fprintf ( 1, '  Set up and execute a job of independent tasks.\n' );
  fprintf ( 1, '  The computation of an integral over [0,1] is divided\n' );
  fprintf ( 1, '  into tasks integrating from A to B.\n' );
  fprintf ( 1, '  The integral estimate is the sum of the task results.\n' );
%
%  Define the job.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Define the job:\n' );

  job = createJob ( ...
    'configuration', 'local', ...
    'FileDependencies', { 'quad_task.m' } );

  task_num = 10;
  n = 100001;
  n_per_task = 1 + ( n - 1 ) / task_num;
%
%  Define the tasks that constitute the job.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Define %d tasks:\n', task_num );

  for task = 1 : task_num
    a = ( task - 1 ) / task_num;
    b =   task       / task_num;
    task_id = createTask ( job, @quad_task, 1, { n, a, b } );
  end
%
%  Submit the job.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Submit the job:\n' );

  submit ( job );
%
%  Wait for results.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Wait for job to finish:\n' );

  timestamp ( )
  wait ( job );
  timestamp ( )
%
%  The following commands will capture and print error messages
%  from the tasks, if something bad occurs.
%
  errmsgs = get ( job.Tasks, {'ErrorMessage'} );
  nonempty = ~cellfun ( @isempty, errmsgs );
  celldisp ( errmsgs(nonempty) );
%
%  Retrieve the output from all tasks.
%
  results = getAllOutputArguments ( job );
%
%  The output is a cell array.
%  We much prefer a MATLAB matrix!
%
  results = cell2mat ( results );
%
%  Sum the results to get the estimate.
%
  quad = sum ( results(:,1) );
  exact = pi;
  err = abs ( quad - exact );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Exact integral = %f\n', exact );
  fprintf ( 1, '  Estimate =       %f\n', quad );
  fprintf ( 1, '  Error =          %e\n', err );
%
%  Clean up by freeing memory associated with the job.
%
  destroy ( job );
%
%  Terminate;
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'QUAD_JOB_LOCAL:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

author: johannesgerer

repo: jburkardt-m


text:
function res = isFeasible(route,model)
    len = length(route);
    if route(1)>model.city||route(len)>model.city
        res = 0; return
    end

    for i = 2:len
        if route(i)>model.city && route(i-1)>model.city  % not use all vehicles
            res = 0; return
        end
    end

    res = 1;
end

author: lzane

repo: VRP-using-SA-with-Matlab


text:
function q = permute(q,varargin)
% overloads permute

q.a = permute(q.a,varargin{:});
q.b = permute(q.b,varargin{:});
q.c = permute(q.c,varargin{:});
q.d = permute(q.d,varargin{:});

author: mtex-toolbox

repo: mtex

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


text:
function varargout = findContour(img, varargin)
%FINDCONTOUR Chain neigbor pixels to form a contour
%
%   Deprecated: replaced by function 'imChainPixels'
%
%   usage:
%   POINTS = findContour(IMG), or
%   [PX, PY] = findContour(IMG)
%   returns a list of points, chaining points visible in the image.
%   IMG shoud be a binary image, containing a single 8-connected loop.
%
%   See Also
%   bwboundaries
%
%   -----
%
%   author: David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 01/11/2003.
%

%   HISTORY
%   06/04/2004 add to graph lib, correct x-y ordering, and add doc.

warning('malimpa:deprecated', ...
    '"findContour" is deprecated, use "imChainPixels" instead');

% set to logical
img = img~=0;

% find points in the image
[pty, ptx] = find(img);

% Initialize iteration : first set the initial point
points = zeros(length(ptx), 2);
x0 = ptx(1);
y0 = pty(1);
points(1, 1) = x0;
points(1, 2) = y0;

% Then set the second point, in the neighbourhood of the first point
vois = findNeighbors(img, [x0 y0]);
x = vois(1, 1);
y = vois(1, 2);
points(2, 1) = x;
points(2, 2) = y;

% For each point, find all neighbours (should be only 2)
% Compare with point previously processed, and choose the other one.
% then shift points references, and loop until all points are processed.
for p=3:length(points)
    vois = findNeighbors(img, [x y]);
    if vois(1,1)==x0 && vois(1,2)==y0
        x0 = x;
        y0 = y;
        x = vois(2, 1);
        y = vois(2, 2);
    else
        x0 = x;
        y0 = y;
        x = vois(1, 1);
        y = vois(1, 2);
    end
    points(p, 1) = x;
    points(p, 2) = y;
end

% format results to match output
if nargout==1
    varargout{1} = points;
elseif nargout==2
    varargout{1} = points(:, 1);
    varargout{2} = points(:, 2);
end

return

function neighList = findNeighbors(img, coord)

xp = coord(1);
yp = coord(2);
neighList = [];
nv = 0;

for x = xp-1:xp+1
    if img(yp-1, x)
        nv = nv+1;
        neighList(nv, 1) = x; %#ok<AGROW>
        neighList(nv, 2) = yp-1;%#ok<AGROW>
    end
    if img(yp+1, x)
        nv = nv+1;
        neighList(nv, 1) = x;%#ok<AGROW>
        neighList(nv, 2) = yp+1;%#ok<AGROW>
    end
end

if img(yp, xp-1)
    nv = nv+1;
    neighList(nv, 1) = xp-1;
    neighList(nv, 2) = yp;
end
if img(yp, xp+1)
    nv = nv+1;
    neighList(nv, 1) = xp+1;
    neighList(nv, 2) = yp;
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
% PREPARE_DATABASE Calculates the features from objects in a source
%
% Usage
%    database = PREPARE_DATABASE(src, feature_fun, options)
%
% Input
%    src (struct): The source specifying the objects.
%    feature_fun (cell): The feature functions applied to each object.
%    options (struct): Options for calculating the features:
%       options.feature_sampling (int): specifies how to sample the feature
%           vectors in time/space (default 1).
%       options.file_normalize (int): The normalization of each file before
%          being given to feature_fun. Can be empty, 1, 2, or Inf (default []).
%       options.parallel (boolean): If true, tries to use the Distributed
%          Computing Toolbox to speed up calculation (default true).
%       Other options are listed in the help for the FEATURE_WRAPPER function.
% Output
%    database (struct): The database of feature vectors.
%
% See also
%    CREATE_SRC, FEATURE_WRAPPER

function db = prepare_database(src,feature_fun,opt)
	if nargin < 3
		opt = struct();
	end

	opt = fill_struct(opt, 'feature_sampling', 1);
	opt = fill_struct(opt, 'file_normalize', []);
	opt = fill_struct(opt, 'parallel', 1);

	features = cell(1,length(src.files));
	obj_ind = cell(1,length(src.files));

	rows = 0;
	cols = 0;
	precision = 'double';

	if opt.parallel
		% parfor loop - note that the contents are the same as the serial loop, but
		% MATLAB doesn't seem to offer an easy way of deduplicating the code.

		% Loop through all the files in the source
		parfor k = 1:length(src.files)
			% Identify all the objects contained in the current file
			file_objects = find([src.objects.ind]==k);

			if length(file_objects) == 0
				% No objects, skip.
				continue;
			end

			tm0 = tic;

			% Load the complete file and normalize as needed.
			x = data_read(src.files{k});

			if ~isempty(opt.file_normalize)
				if opt.file_normalize == 1
					x = x/sum(abs(x(:)));
				elseif opt.file_normalize == 2
					x = x/sqrt(sum(abs(x(:)).^2));
				elseif opt.file_normalize == Inf
					x = x/max(abs(x(:)));
				end
			end

			% Due to parallelization constraints, each file stores its feature vectors
			% in one index of the features cell array. The corresponding object indices
			% are stored in the obj_ind cell array.
			features{k} = cell(1,length(file_objects));
			obj_ind{k} = file_objects;

			% Apply all the feature functions to the objects contained in x. The output
			% buf is a PxNxK vector, where P corresponds to number of feature coeffi-
			% cients, N corresponds to the number of features - the number of "windows"
			% - in each object, and K corresponds to the number of objects.
			buf = apply_features(x,src.objects(file_objects),feature_fun,opt);

			% Subsample among the features as needed.
			buf = buf(:,1:opt.feature_sampling:end,:);

			% Separate each of the objects into its own array element.
			for l = 1:length(file_objects)
				features{k}{l} = buf(:,:,l);
			end
			fprintf('.');
		end
	else
		time_start = clock;
		for k = 1:length(src.files)
			file_objects = find([src.objects.ind]==k);

			if length(file_objects) == 0
				continue;
			end

			tm0 = tic;
			x = data_read(src.files{k});

			if ~isempty(opt.file_normalize)
				if opt.file_normalize == 1
					x = x/sum(abs(x(:)));
				elseif opt.file_normalize == 2
					x = x/sqrt(sum(abs(x(:)).^2));
				elseif opt.file_normalize == Inf
					x = x/max(abs(x(:)));
				end
			end

			features{k} = cell(1,length(file_objects));
			obj_ind{k} = file_objects;

			buf = apply_features(x,src.objects(file_objects),feature_fun,opt);

			buf = buf(:,1:opt.feature_sampling:end,:);

			for l = 1:length(file_objects)
				features{k}{l} = buf(:,:,l);
			end
			time_elapsed = etime(clock, time_start);
			estimated_time_left = time_elapsed * (length(src.files)-k) / k;
			fprintf('.');
		end
	end
	fprintf('\n');

	% Identify the files for which objects have been computed.
	nonempty = find(~cellfun(@isempty,features));

	% Determine the size of the features array. Each row corresponds to a feature
	% coefficient, each column corresponds to a feature vector.
	rows = size(features{nonempty(1)}{1},1);
	cols = sum(cellfun(@(x)(sum(cellfun(@(x)(size(x,2)),x))),features(nonempty)));
	precision = class(features{nonempty(1)}{1});

	db.src = src;

	db.features = zeros(rows,cols,precision);

	db.indices = cell(1,length(src.objects));

	% Fill in the features array using the pre-calculated features.
	r = 1;
	for k = 1:length(features)
		for l = 1:length(features{k})
			ind = r:r+size(features{k}{l},2)-1;
			db.features(:,ind) = features{k}{l};
			db.indices{obj_ind{k}(l)} = ind;
			r = r+length(ind);
		end
	end
end

function out = apply_features(x,objects,feature_fun,opt)
	% For each feature function, apply it to the objects in x and concatenate
	% the results.

	if ~iscell(feature_fun)
		feature_fun = {feature_fun};
	end

	out = {};
	for k = 1:length(feature_fun)
		if nargin(feature_fun{k}) == 2
			% If only two inputs are given, the feature function will handle the cut
			% ting up of x into objects itself.
			out{k,1} = feature_fun{k}(x,objects);
		elseif nargin(feature_fun{k}) == 1
			% If only one input is given, the feature function relies on us to cut up
			% the input x into objects and feed them to it.
			out{k,1} = feature_wrapper(x,objects,feature_fun{k},opt);
		else
			error('Invalid number of inputs for feature function!')
		end
	end

	if length(feature_fun) == 1
		out = out{1};
	else
		% Concatenate along the first dimension - the feature dimension - of each
		% element.
		out = cellfun(@(x)(permute(x,[2 1 3])),out,'UniformOutput',false);
		out = permute([out{:}],[2 1 3]);
	end
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function my_net_run_batch(net_config, work_info_batch)

    gen_data_info_groups_cache(net_config, work_info_batch);
    group_idx=net_config.ref.root_group_idx;

    data_info=init_one_data_info_from_cache(group_idx, work_info_batch);
    input_info=work_info_batch.ref.ds_info.net_input_info;
    data_info.ref.output_info_layers{1}=input_info;
    work_info_batch.ref.ds_info.net_input_info=[];

    do_run_one_pass(net_config, work_info_batch, group_idx);

    work_info_batch.ref.data_info_groups=[];

end

function do_run_one_pass(net_config, work_info_batch, group_idx)

    if work_info_batch.ref.net_run_verbose
        fprintf('---- runing one pass, root_group_idx:%d\n', group_idx);
    end

    do_run_one_group(net_config, work_info_batch, group_idx, 'forward');

    data_info=work_info_batch.ref.data_info_groups{group_idx};
    last_layer_output_info=data_info.ref.output_info_layers{end};
    assert(check_valid_net_output(last_layer_output_info));

    if work_info_batch.ref.run_backward

        init_dzdx=single(1);

        last_layer_output_info.dzdx=init_dzdx;
        last_layer_output_info.bp_finished=true;
        data_info.ref.output_info_layers{end}=last_layer_output_info;

        do_run_one_group(net_config, work_info_batch, group_idx, 'backward');
    end

end

function do_run_one_group(net_config, work_info_batch, group_idx, group_run_type)

group_info=net_config.ref.group_infos{group_idx};
work_info_batch.ref.net_run_current_group_idx=group_idx;

net_run_current_group_idx_linked=[];
if ~isempty(group_info.net_info)
    if isfield(group_info.net_info.ref, 'linked_group_idx')
        net_run_current_group_idx_linked=group_info.net_info.ref.linked_group_idx;
    end
end
work_info_batch.ref.net_run_current_group_idx_linked=net_run_current_group_idx_linked;

if work_info_batch.ref.net_run_verbose
    fprintf('----start group_idx:%d, name:%s, run_type:%s\n', ...
        group_info.group_idx, group_info.name, group_run_type);
end

if strcmp(group_run_type, 'forward')
    do_forward_one_group( net_config, work_info_batch, group_idx);
end

if strcmp(group_run_type, 'backward')
    do_backward_one_group( net_config, work_info_batch, group_idx);
end

if work_info_batch.ref.net_run_verbose
    fprintf('----finish group_idx:%d, name:%s, run_type:%s\n', ...
        group_info.group_idx, group_info.name, group_run_type);
end

end

function do_forward_group_child_chain( work_info_batch, data_info, group_info)

        net_config=work_info_batch.ref.net_config;

        child_group_idxes=group_info.child_group_idxes;

        child_num=length(child_group_idxes);
        child_valid_flags=data_info.ref.child_valid_flags;
        assert(any(child_valid_flags));

        input_info=data_info.ref.output_info_layers{1};
        assert(check_valid_net_output(input_info));
        need_bp=false;

        for g_idx=1:child_num

            child_valid=child_valid_flags(g_idx);
            if ~child_valid
                break;
            end

            child_group_idx=child_group_idxes(g_idx);

            [input_info, need_bp]=do_forward_one_child( net_config, ...
                work_info_batch, child_group_idx, input_info);

            if ~check_valid_net_output(input_info)
                break;
            end
        end

        data_info.ref.output_info_layers{2}=input_info;
        data_info.ref.need_bp=need_bp;

end

function do_forward_group_child_parallel( work_info_batch, data_info, group_info)

    net_config=work_info_batch.ref.net_config;

    child_group_idxes=group_info.child_group_idxes;
    child_num=length(child_group_idxes);
    child_valid_flags=data_info.ref.child_valid_flags;
    assert(any(child_valid_flags));

    group_input_info=data_info.ref.output_info_layers{1};

    assert(check_valid_net_output(group_input_info));

    assert(group_input_info.is_group_data);

    input_info_child_groups=group_input_info.data_child_groups;
    input_child_valid_flags=group_input_info.child_valid_flags;
    assert(length(input_info_child_groups)==child_num);
    assert(length(input_child_valid_flags)==child_num);

    output_child_groups=cell(child_num, 1);
    output_valid_flags=false(child_num,1);
    need_bp_flags=false(child_num,1);

    for g_idx=1:child_num

        child_valid=child_valid_flags(g_idx);
        if ~child_valid
            continue;
        end

        input_valid=input_child_valid_flags(g_idx);
        if ~input_valid
            continue;
        end
        child_input_info=input_info_child_groups{g_idx};

        child_group_idx=child_group_idxes(g_idx);

        [child_output_info, child_need_bp]=do_forward_one_child(net_config, ...
            work_info_batch, child_group_idx, child_input_info);
        output_child_groups{g_idx}=child_output_info;

        need_bp_flags(g_idx)=child_need_bp;

        output_valid_flags(g_idx)=check_valid_net_output(child_output_info);
    end

    assert(any(output_valid_flags));

    group_output_info=[];
    group_output_info.data_child_groups=output_child_groups;
    group_output_info.child_valid_flags=output_valid_flags;
    group_output_info.is_group_data=true;
    group_output_info.forward_finished=any(output_valid_flags);

    data_info.ref.output_info_layers{2}=group_output_info;
    data_info.ref.need_bp=any(need_bp_flags);
    data_info.ref.need_bp_child_groups=need_bp_flags;

end

function do_forward_one_group( net_config, work_info_batch, group_idx)

group_info=net_config.ref.group_infos{group_idx};
prediction_layer_idxes=group_info.prediction_layer_idxes;

if ~isempty(group_info.forward_begin_fn)
    group_info=group_info.forward_begin_fn(work_info_batch, group_info);
end

if ~group_info.skip_forward

    child_group_idxes=group_info.child_group_idxes;
    data_info=work_info_batch.ref.data_info_groups{group_idx};

    if isempty(child_group_idxes)

        net_run_config=work_info_batch.ref.net_run_config;

        extra_output_layer_idxes=[];
        if ~isempty(prediction_layer_idxes)
            extra_output_layer_idxes=cat(1, extra_output_layer_idxes, prediction_layer_idxes);
            extra_output_layer_idxes=unique(extra_output_layer_idxes);
        end

        work_info_batch.ref.current_group_data_info=data_info;

        forward_time=tic;
        extra_data_info=my_net_forward(retrieve_net_info(group_info, net_config), ...
            work_info_batch, data_info, net_run_config, extra_output_layer_idxes);
        forward_time=toc(forward_time);
        work_info_batch.ref.forward_time=work_info_batch.ref.forward_time+forward_time;

        data_info.ref.extra_input_infos=[];

        work_info_batch.ref.current_group_data_info=[];

    else

        switch group_info.child_relation
            case 'chain'
                do_forward_group_child_chain( work_info_batch, data_info, group_info);
            case 'parallel'
                do_forward_group_child_parallel( work_info_batch, data_info, group_info);
            otherwise
                error('not support!');
        end

    end
end

if ~isempty(prediction_layer_idxes)

    assert(isempty(child_group_idxes));

    prediction_info=work_info_batch.ref.prediction_info_groups{group_idx};
    keep_output_infos=extra_data_info.output_layers(group_info.prediction_layer_idxes);

    prediction_info.output_infos=keep_output_infos;
    work_info_batch.ref.prediction_info_groups{group_idx}=prediction_info;
end

if ~isempty(group_info.forward_finish_fn)
    group_info=group_info.forward_finish_fn(work_info_batch, group_info);
end

if ~isempty(group_info.forward_evaluate_fn)
    forward_eva_time=tic;
    group_info.forward_evaluate_fn( work_info_batch, group_idx);
    forward_eva_time=toc(forward_eva_time);
    work_info_batch.ref.eva_time=work_info_batch.ref.eva_time+forward_eva_time;
end

end

function [output_info, child_need_bp]=do_forward_one_child(net_config,...
    work_info_batch, child_group_idx, input_info)

    child_data_info=init_one_data_info_from_cache(child_group_idx, work_info_batch);
    child_data_info.ref.output_info_layers{1}=input_info;

    do_run_one_group(net_config, work_info_batch, child_group_idx, 'forward');
    output_info=child_data_info.ref.output_info_layers{end};

    child_data_info.ref.output_info_layers{end}=[];
    child_need_bp=child_data_info.ref.need_bp;

    if ~child_need_bp
        child_data_info.ref.output_info_layers=[];
    end

end

function output_info_layers=gen_net_output_info_basic(net_info)

if isempty(net_info)
    n=1;
else
    n = numel(net_info.ref.layers) ;
end

sample_output_info=[];
sample_output_info.x=[];
sample_output_info.aux=[];
sample_output_info.dzdx=[];
sample_output_info.dzdw=[];
sample_output_info.is_group_data=false;
sample_output_info.bp_finished=false;
sample_output_info.forward_finished=false;

output_info_layers=repmat({sample_output_info}, n+1, 1);

end

function data_info=init_one_data_info_from_cache(group_idx, work_info_batch)

data_info=work_info_batch.ref.data_info_groups{group_idx};
data_info.ref=work_info_batch.ref.sample_data_info_groups{group_idx};

tmp_data_info=work_info_batch.ref.data_info_groups_extra{group_idx};
data_info.ref.extra_input_infos=tmp_data_info.extra_input_infos;

end

function data_info=init_one_data_info(group_info, net_config)

    data_info=[];
    data_info.output_info_layers=[];
    data_info.forward_time=0;
    data_info.backward_time=0;
    data_info.need_bp=false;
    data_info.extra_input_infos=cell(0);

    child_num=length(group_info.child_group_idxes);
    data_info.child_valid_flags=true(child_num, 1);

    data_info.output_info_layers=gen_net_output_info_basic(...
        retrieve_net_info(group_info, net_config));

end

function gen_data_info_groups_cache(net_config, work_info_batch)

    group_num=length(net_config.ref.group_infos);
    work_info=work_info_batch.ref.work_info;
    if ~isfield(work_info.ref.tmp_cache, 'sample_data_info_groups')
        data_info_groups=cell(group_num, 1);
        data_info_groups_ref=cell(group_num, 1);
        for g_idx=1:group_num
            tmp_data_info=[];
            make_ref_obj(tmp_data_info);
            data_info_groups_ref{g_idx}=tmp_data_info;
            data_info=init_one_data_info(net_config.ref.group_infos{g_idx}, net_config);
            data_info_groups{g_idx}=data_info;
        end
        work_info.ref.tmp_cache.sample_data_info_groups=data_info_groups;
        work_info.ref.tmp_cache.sample_data_info_groups_ref=data_info_groups_ref;

        data_info_groups_extra=cell(group_num, 1);
        for g_idx=1:group_num
            tmp_data_info_extra=[];
            tmp_data_info_extra.extra_input_infos=cell(0);
            data_info_groups_extra{g_idx}=tmp_data_info_extra;
        end
        work_info.ref.tmp_cache.sample_data_info_groups_extra=data_info_groups_extra;

    end

    work_info_batch.ref.sample_data_info_groups=work_info.ref.tmp_cache.sample_data_info_groups;
    work_info_batch.ref.data_info_groups=work_info.ref.tmp_cache.sample_data_info_groups_ref;
    work_info_batch.ref.prediction_info_groups=cell(group_num, 1);

    work_info_batch.ref.data_info_groups_extra=work_info.ref.tmp_cache.sample_data_info_groups_extra;

end

function do_backward_group_child_chain(net_config,work_info_batch, data_info, group_info)

    child_group_idxes=group_info.child_group_idxes;
    child_num=length(child_group_idxes);
    child_valid_flags=data_info.ref.child_valid_flags;
    assert(any(child_valid_flags));

        start_g_idx=find(~child_valid_flags, 1);
        if isempty(start_g_idx)
            start_g_idx=child_num;
        end

        init_output_info=data_info.ref.output_info_layers{end};

        for g_idx=start_g_idx:-1:1

            child_group_idx=child_group_idxes(g_idx);
            init_output_info=do_backward_child_group( net_config,...
                work_info_batch, child_group_idx, init_output_info);

            if ~init_output_info.bp_finished
                break;
            end
        end

        data_info.ref.output_info_layers{1}=init_output_info;

end

function do_backward_group_child_parallel(net_config,work_info_batch,data_info, group_info)

    child_group_idxes=group_info.child_group_idxes;
    child_num=length(child_group_idxes);
    child_valid_flags=data_info.ref.child_valid_flags;
    assert(any(child_valid_flags));

        init_group_output_info=data_info.ref.output_info_layers{end};

        assert(init_group_output_info.is_group_data);
        output_info_child_groups=init_group_output_info.data_child_groups;
        output_child_valid_flags=init_group_output_info.child_valid_flags;
        init_bp_child_valid_flags=init_group_output_info.bp_child_valid_flags;

        assert(length(output_info_child_groups)==child_num);
        assert(length(output_child_valid_flags)==child_num);

        first_output_info=data_info.ref.output_info_layers{1};
        assert(first_output_info.is_group_data);

        bp_child_valid_flags=false(child_num, 1);
        for g_idx=1:child_num

            child_valid=output_child_valid_flags(g_idx);
            if ~child_valid
                continue;
            end

            init_bp_finished=init_bp_child_valid_flags(g_idx);
            if ~init_bp_finished
                continue;
            end

            child_group_idx=child_group_idxes(g_idx);
            child_init_output_info=output_info_child_groups{g_idx};

            if ~check_valid_net_output(child_init_output_info)
                continue;
            end

            child_first_output_info=do_backward_child_group( net_config,...
                work_info_batch, child_group_idx, child_init_output_info);

            first_output_info.data_child_groups{g_idx}=child_first_output_info;

            bp_child_valid_flags(g_idx)=child_first_output_info.bp_finished;

        end

        first_output_info.bp_child_valid_flags=bp_child_valid_flags;
        first_output_info.bp_finished=any(bp_child_valid_flags);

        data_info.ref.output_info_layers{1}=first_output_info;
end

function do_backward_one_group( net_config, work_info_batch, group_idx)

data_info=work_info_batch.ref.data_info_groups{group_idx};
if ~data_info.ref.need_bp
    return;
end

group_info=net_config.ref.group_infos{group_idx};

if ~group_info.skip_backward

    child_group_idxes=group_info.child_group_idxes;

    if isempty(child_group_idxes)

        net_run_config=work_info_batch.ref.net_run_config;

        backward_time=tic;
        do_bp_current_net=check_do_bp_current_net(work_info_batch, net_run_config, ...
            retrieve_net_info(group_info, net_config));
        if do_bp_current_net

            work_info_batch.ref.current_group_data_info=data_info;

            my_net_backward(retrieve_net_info(group_info, net_config), ...
                work_info_batch, data_info, net_run_config);
            work_info_batch.ref.leaf_group_bp_flags(group_idx)=true;

            work_info_batch.ref.current_group_data_info=[];

        end

        backward_time=toc(backward_time);
        work_info_batch.ref.backward_time=work_info_batch.ref.backward_time+backward_time;

    else

        switch group_info.child_relation
            case 'chain'
                do_backward_group_child_chain(net_config, work_info_batch, data_info, group_info);
            case 'parallel'
                do_backward_group_child_parallel(net_config, work_info_batch, data_info, group_info);
            otherwise
                error('not support!');
        end

    end
end

data_info.ref.output_info_layers=data_info.ref.output_info_layers(1);

end

function first_output_info=do_backward_child_group( ...
    net_config, work_info_batch, child_group_idx, init_output_info)

    child_data_info=work_info_batch.ref.data_info_groups{child_group_idx};

    if ~child_data_info.ref.need_bp
        first_output_info=[];
        first_output_info.bp_finished=false;
        return;
    end

    if ~check_valid_net_output(init_output_info)
        first_output_info=[];
        first_output_info.bp_finished=false;
        return;
    end

    child_data_info.ref.output_info_layers{end}=init_output_info;
        do_run_one_group(net_config, work_info_batch, child_group_idx, 'backward');

    first_output_info=child_data_info.ref.output_info_layers{1};
    if isempty(first_output_info)
        first_output_info.bp_finished=false;
    end

    child_data_info.ref=[];

end

function one_net_info=retrieve_net_info(group_info, net_config)

one_net_info=group_info.net_info;
if isempty(one_net_info)
    return;
end
if isempty(one_net_info.ref.layers)
    linked_group_info=net_config.ref.group_infos{one_net_info.ref.linked_group_idx};
    one_net_info=linked_group_info.net_info;
end

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
% Flush out the MATLAB.
close all;
clc;
clear all;

% Read the desired image file.
ImageData=imread('building.pnm');

% Add noise to oriiginal image.
ImageData = imnoise(ImageData,'Salt & Pepper', 0.04);

% Display the original image.
figure,imshow(ImageData);
title(' Original Image with noise: ');

% Define constant.
Constant=0.3;

% Create the matrix.
[X,Y]=size(ImageData);

        for a = 1:X

            for b = 1:Y
                m=double(ImageData(a,b));

                OutputImage(a,b)=Constant.*log10(1+m);
            end

        end

% Display the output image.
figure, imshow(OutputImage);
title(' Final Image: ');

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function result = sphere_monomial_int_nd ( n, r, e )

%*****************************************************************************80
%
%% SPHERE_MONOMIAL_INT_ND integrates a monomial on a sphere in ND.
%
%  Discussion:
%
%    The sphere may have nonunit radius, but it must be centered at 0.
%
%    The integration region is
%
%      sum ( x(i)**2 ) = R**2.
%
%    The monomial is F(X) = X(1)**E(1) * X(2)**E(2) * ... * X(N)**E(N).
%
%    This routine is useful for testing the accuracy of quadrature
%    rules on the sphere.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 January 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Philip Davis, Philip Rabinowitz,
%    Methods of Numerical Integration,
%    Second Edition,
%    Academic Press, 1984, page 263.
%
%  Parameters:
%
%    Input, integer N, the dimension of the space.
%
%    Input, real R, the radius of the sphere.
%
%    Input, integer E(N), the exponents of X, Y and Z in the monomial.
%    Each exponent must be nonnegative.
%
%    Output, real RESULT, the integral.
%
  if ( any ( e(1:n) < 0 ) )
    result = -r8_huge ( result );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SPHERE_MONOMIAL_INT_ND - Fatal error!\n' );
    fprintf ( 1, '  All exponents must be nonnegative.\n' );
    error ( 'SPHERE_MONOMIAL_INT_ND - Fatal error!' );
  end

  if ( all ( e(1:n) == 0 ) )

    result = 2.0 * sqrt ( pi^n ) / gamma ( 0.5 * n );

  elseif ( any ( mod ( e(1:n), 2 ) == 1 ) )

    result = 0.0;

  else

    result = 2.0;

    for i = 1 : n
      result = result * gamma ( 0.5 * ( e(i) + 1 ) );
    end

    result = result / gamma ( 0.5 * ( sum ( e(1:n) + 1 ) ) );

  end

  result = result * r^( sum ( e(1:n) ) + 2 );

  return
end


text:
function x = creation_null_right ( m, n )

%*****************************************************************************80
%
%% CREATION_NULL_RIGHT returns a right null vector of the CREATION matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the order of A.
%
%    Output, real X(N,1), the null vector.
%
  x = zeros ( n, 1 );

  x(n,1) = 1.0;

  return
end


text:
function [N,errL2,errH1,erruIuh] = cubePoissoncvtmeshnew
% CUBEPOISSON solves Poisson equation in a cube.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

close all; clear all;
%% Parameters
maxIt = 4; N = zeros(maxIt,1);
errL2 = zeros(maxIt,1); errH1 = zeros(maxIt,1); erruIuh = zeros(maxIt,1);

%% Get the data of the pde
pde = sincosdata3;
% pde = polydata3;

%% Finite Element Method
for k = 1:maxIt
    [node,elem] = cvtuniformmesh(2^k);
    % solve the equation
    option.solver = 'direct';
    [u,Du,A,eqn] = Poisson3(node,elem,pde,[],[],option);
    N(k) = size(node,1);
    % compute error
    errL2(k) = getL2error3(node,elem,pde.exactu,u);
    errH1(k) = getH1error3(node,elem,pde.Du,Du);
    uI = pde.exactu(node);  % nodal interpolation
    erruIuh(k) = sqrt((u-uI)'*eqn.AD*(u-uI));
end

%% Plot convergence rates
figure(2);
r1 = showrate(N,errH1,2,'-*');
hold on;
r2 = showrate(N,errL2,2,'k-+');
r3 = showrate(N,erruIuh,2,'m-+');
legend('||Du-Du_h||',['N^{' num2str(r1) '}'], ...
       '||u-u_h||',['N^{' num2str(r2) '}'], ...
       '||DuI-Du_h||',['N^{' num2str(r3) '}'], 'LOCATION','Best');
%%
% The error in H1 and L2 norm converges at optimal rate. But no
% superconvergence on criss-cross grids.
end

repo: Aero-matlab


text:
function [TorF, vstr, rdate] = have_feature_pdipmopf()
%HAVE_FEATURE_PDIPMOPF  Detect availability/version info for PDIPMOPF
%
%   Feature detection function implementing 'pdipmopf' tag for HAVE_FEATURE
%   to detect availability/version of PDIPMOPF, a primal-dual interior point
%   method optimal power flow (OPF) solver included in TSPOPF.
%   (https://www.pserc.cornell.edu/tspopf)
%
%   See also HAVE_FEATURE, PDIPMOPF.

%   MATPOWER
%   Copyright (c) 2004-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

if have_feature('matlab') && exist('pdipmopf', 'file') == 3;
    TorF = 1;
    v = pdipmopfver('all');
    vstr = v.Version;
    rdate = v.Date;
else
    TorF = 0;
    vstr = '';
    rdate = '';
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function simon(varargin)
%SIMON Do what Simon Says...
%   SIMON starts the standard game.
%   SIMON(N) starts the game with N keys.
%
%   Instructions:
%   - Press the start button (the green triangle pointing to the right) to
%     start the game.
%   - Look and listen to what Simon says and repeat it.
%   - Press the stop button (the red square) to stop the game.
%   - Press the up and down buttons (the yellow triangles pointing up and
%     down) to change the number of buttons.
%   - The numbers in the middle give the current score, the high score and
%     the numbers of buttons.
%
%   Example:
%     simon    start the standard game with four buttons.
%     simon(6) start the game with six buttons.

%   Eduard van der Zwan, 21-Jan-2007

%handle input
if nargin==0
    numberOfButtons=4;
elseif nargin==1
    numberOfButtons=varargin{1};
    if ~isnumeric(numberOfButtons)||rem(numberOfButtons,1)~=0||numberOfButtons<=1
        error('Input must be an integer larger then 1.')
    end
else
    error('Too many input arguments.')
end

%load high score file
highScoreFile=[fileparts(which(mfilename)) filesep 'simonHighScore.mat'];
if exist(highScoreFile,'file')
    load(highScoreFile)
    if length(highScore)<numberOfButtons
        %we started with a high number of buttons
        highScore(numberOfButtons)=0;
    end
else
    highScore=zeros(numberOfButtons,1);
end
%create figure
h.figure=figure(...
    'Name','Simon',...
    'NumberTitle','off',....
    'Color','k',...
    'Toolbar','none',...
    'MenuBar','none',...
    'CloseRequestFcn',@closeRequest);

%user data
g=struct(...
    'buttonTime',0.2,...seconds
    'intervalTime',0.3,...seconds
    'sampleFrequency',10000,...Hz
    'dimRatio',0.5,...ratio between lid and dimmed button color
    'innerRadius',1.2,...of the big buttons
    'outerRadius',3,...of the big buttons
    'buttonStepRatio',0.1,...ratio between gap and button size
    'startButtonColor',[0 1 0],...green
    'stopButtonColor',[1 0 0],...red
    'upDownButtonColor',[1 1 0],...yellow
    'buttonRadius',0.2,...of the small button
    'oldPosition',get(h.figure,'Position'),...
    'numberOfButtons',numberOfButtons,...
    'highScore',highScore,...
    'highScoreFile',highScoreFile,...
    'userSequence',[],...
    'simonSequence',[],...
    'playUserSequence',false);

set(h.figure,'UserData',g);

%create axes
h.axes=axes(...
    'DataAspectRatio',[1 1 1],...
    'NextPlot','add',...
    'Visible','off',...
    'Parent',h.figure);

%create small buttons
%location of buttuns relative to origin
xButton=[-2 -2 2 2]*g.buttonRadius;
yButton=[2 -2 2 -2]*g.buttonRadius;

%create start button
angle=linspace(0,2*pi,4);
x=xButton(1)+cos(angle)*g.buttonRadius;
y=yButton(1)+sin(angle)*g.buttonRadius;
h.startButton=patch(...
    x,y,g.startButtonColor,...
    'EdgeColor','none',...
    'Parent',h.axes,...
    'BusyAction','cancel',...
    'Interruptible','off');

%create stop button
angle=linspace(0.25*pi,2.25*pi,5);
x=xButton(2)+cos(angle)*g.buttonRadius;
y=yButton(2)+sin(angle)*g.buttonRadius;
h.stopButton=patch(...
    x,y,g.stopButtonColor*g.dimRatio,...
    'EdgeColor','none',...
    'Parent',h.axes,...
    'BusyAction','cancel',...
    'Interruptible','off');

%create up button
angle=linspace(0.5*pi,2.5*pi,4);
x=xButton(3)+cos(angle)*g.buttonRadius;
y=yButton(3)+sin(angle)*g.buttonRadius;
h.upButton=patch(...
    x,y,g.upDownButtonColor,...
    'EdgeColor','none',...
    'Parent',h.axes,...
    'BusyAction','cancel',...
    'Interruptible','off');

%create down button
angle=linspace(1.5*pi,3.5*pi,4);
x=xButton(4)+cos(angle)*g.buttonRadius;
y=yButton(4)+sin(angle)*g.buttonRadius;
h.downButton=patch(...
    x,y,g.upDownButtonColor,...
    'EdgeColor','none',...
    'Parent',h.axes,...
    'BusyAction','cancel',...
    'Interruptible','off');

%create score text
h.score=text(...
    'Parent',h.axes,...
    'HorizontalAlignment','Center',...
    'VerticalAlignment','Middle',...
    'FontSize',18,...
    'FontWeight','Bold',...
    'Color','w');
setScore(h)

%create buttons in circle
createButtons(h)

function createButtons(h)
%create buttons in the circle
g=get(h.figure,'UserData');

%determine colors and sounds
%about sounds:
%   the 'A' has a frequency of 440 Hz. All other frequncies of the piano
%   keys can be be determined with:
%       keyFrequency=440*2^(keyNumber/12)
%   in which keyNumber is the number of the key relative to the 'A' key;
%   with 'A' being 0, left of that key minus numbers, and right of that key
%   plus numbers
if g.numberOfButtons==4
    %original colors and sounds
    buttonCollor=[1 0 0;0 1 0;1 1 0;0 0 1]*g.dimRatio;
    %red = [1 0 0], green = [0 1 0], yellow = [1 1 0] and blue = [0 0 1]
    buttonFrequency=440*2.^([0 12 10 5]/12);
    %A = 0, a = 12, G = 10 and D = 5
else
    buttonCollor=hsv(g.numberOfButtons)*g.dimRatio;
    buttonFrequency=440*2.^((0:g.numberOfButtons-1)/12);
end

buttonWidthAngle=2*pi/g.numberOfButtons;
outerAngleStep=buttonWidthAngle*g.buttonStepRatio;
innerAngleStep=outerAngleStep*g.outerRadius/g.innerRadius;

h.button=[];
for buttonNumber=1:g.numberOfButtons
    %user data for button
    signalVector=sin(linspace(0,buttonFrequency(buttonNumber)*g.buttonTime*2*pi,g.sampleFrequency*g.buttonTime));
    b=struct(...
        'signalVector',signalVector,...
        'buttonNumber',buttonNumber);
    %draw button
    buttonAngle=(buttonNumber-1)*buttonWidthAngle;
    innerArc=linspace(buttonAngle,buttonAngle+buttonWidthAngle-innerAngleStep,20)+innerAngleStep/2;
    outerArc=linspace(buttonAngle,buttonAngle+buttonWidthAngle-outerAngleStep,40)+outerAngleStep/2;
    x=[g.innerRadius*cos(innerArc) g.outerRadius*cos(fliplr(outerArc))];
    y=[g.innerRadius*sin(innerArc) g.outerRadius*sin(fliplr(outerArc))];
    h.button(buttonNumber)=patch(...
        x,y,buttonCollor(buttonNumber,:),...
        'EdgeColor','none',...
        'Parent',h.axes,...
        'UserData',b,...
        'BusyAction','cancel',...
        'Interruptible','off');
end

%set object functions including the handles to above buttons
for buttonNumber=1:g.numberOfButtons
    set(h.button(buttonNumber),'ButtonDownFcn',{@buttonDown,h})
end
set(h.startButton,'ButtonDownFcn',{@startButtonDown,h});
set(h.stopButton,'ButtonDownFcn',{@stopButtonDown,h});
set(h.upButton,'ButtonDownFcn',{@upButtonDown,h});
set(h.downButton,'ButtonDownFcn',{@downButtonDown,h});
set(h.figure,'ResizeFcn',{@resize,h})

function startButtonDown(hStartButton,eventdata,h)
%start button is pressed
g=get(h.figure,'UserData');
startButtonColor=get(hStartButton,'FaceColor');
if isequal(startButtonColor,g.startButtonColor)
    %make buttons inactive
    set(hStartButton,'FaceColor',g.startButtonColor*g.dimRatio)
    set(h.upButton,'FaceColor',g.upDownButtonColor*g.dimRatio)
    set(h.downButton,'FaceColor',g.upDownButtonColor*g.dimRatio)
    %make stop button active
    set(h.stopButton,'FaceColor',g.stopButtonColor)
    playSimonSequence(h)
end

function stopButtonDown(hStopButton,eventdata,h)
%stop button is pressed
g=get(h.figure,'UserData');
stopButtonColor=get(hStopButton,'FaceColor');
if isequal(stopButtonColor,g.stopButtonColor)
    %make stop button inactive
    set(hStopButton,'FaceColor',g.stopButtonColor*g.dimRatio)
    %make buttons active
    set(h.startButton,'FaceColor',g.startButtonColor)
    set(h.upButton,'FaceColor',g.upDownButtonColor)
    set(h.downButton,'FaceColor',g.upDownButtonColor)
    %reset data
    g.userSequence=[];
    g.simonSequence=[];
    g.playUserSequence=false;
    set(h.figure,'UserData',g)
end

function upButtonDown(hUpButton,eventdata,h)
%up button is pressed
g=get(h.figure,'UserData');
upButtonColor=get(hUpButton,'FaceColor');
if isequal(upButtonColor,g.upDownButtonColor)
    g.numberOfButtons=g.numberOfButtons+1;
    %expand high score
    if length(g.highScore)<g.numberOfButtons
        g.highScore(g.numberOfButtons)=0;
    end
    set(h.figure,'UserData',g)
    setScore(h)
    delete(h.button)
    %create buttons in circle
    createButtons(h)
end

function downButtonDown(hDownButton,eventdata,h)
%down button is pressed
g=get(h.figure,'UserData');
downButtonColor=get(hDownButton,'FaceColor');
if isequal(downButtonColor,g.upDownButtonColor) && g.numberOfButtons>2
    g.numberOfButtons=g.numberOfButtons-1;
    set(h.figure,'UserData',g)
    setScore(h)
    delete(h.button)
    %create buttons in circle
    createButtons(h)
end

function buttonDown(hButton,eventdata,h)
%a button in the circle is pressed
g=get(h.figure,'UserData');
b=get(hButton,'UserData');

dimColor=get(hButton,'FaceColor');
litColor=dimColor/g.dimRatio;
%lid button and play sound
set(hButton,'FaceColor',litColor)
sound(b.signalVector,g.sampleFrequency)
pause(g.buttonTime)
set(hButton,'FaceColor',dimColor)
if g.playUserSequence
    %button is pressed by user
    g.userSequence(end+1)=b.buttonNumber;
    set(h.figure,'UserData',g);
    if isequal(g.userSequence,g.simonSequence)
        %user repeated Simon correctly
        playSimonSequence(h)
    elseif ~isequal(g.userSequence,g.simonSequence(1:length(g.userSequence)))
        %user made mistake
        stopButtonDown(h.stopButton,[],h)
        %play rest of Simons sequence
        restSequence=g.simonSequence(length(g.userSequence):end);
        for counter=1:length(restSequence)
            pause(g.intervalTime)
            hButton=h.button(restSequence(counter));
            buttonDown(hButton,[],h);
        end
    end
end

function closeRequest(hFigure,eventdata)
%figure is closed
g=get(hFigure,'UserData');
%save high score
highScore=g.highScore;
save(g.highScoreFile,'highScore')
%and delete figure
%error if a sequence is played
delete(hFigure)

function resize(hFigure,eventdata,h)
%figure is resized
g=get(hFigure,'UserData');
%change font of score text
newPosition=get(hFigure,'Position');
newFontSize=get(h.score,'FontSize')*min(newPosition(3:4))/min(g.oldPosition(3:4));
set(h.score,'FontSize',newFontSize)
g.oldPosition=newPosition;
set(hFigure,'UserData',g);

function playSimonSequence(h)
%Simon playes his sequence
g=get(h.figure,'UserData');
%change high score
if length(g.simonSequence)>g.highScore(g.numberOfButtons)
    g.highScore(g.numberOfButtons)=length(g.simonSequence);
end
set(h.figure,'UserData',g)
setScore(h)
%expand Simons sequence
g.simonSequence(end+1)=ceil(rand*g.numberOfButtons);
g.userSequence=[];
g.playUserSequence=false;
set(h.figure,'UserData',g)
%play sequence
for counter=1:length(g.simonSequence)
    pause(g.intervalTime)
    hButton=h.button(g.simonSequence(counter));
    buttonDown(hButton,[],h);
end
g.playUserSequence=true;
set(h.figure,'UserData',g);

function setScore(h)
%reset score text
g=get(h.figure,'UserData');
score=int2str(length(g.simonSequence));
highScore=int2str(g.highScore(g.numberOfButtons));
numberOfButtons=int2str(g.numberOfButtons);
set(h.score,'String',[score ' ' highScore ' ' numberOfButtons])

author: Sable

repo: mcbench-benchmarks


text:
function [rho,u,p,fl,Ml_eq] = ModSBBGK()
%function [rho_our p_out t_out]=M_codeA(rho_in,p_in,t_in)
 %switch method
 %   case{1} % TVD 0(h^2)
        % Using discrete ordinate method (discrete and constant velocity
        % values in phase-space domain)
        a = v(:,1);
        % Load initial condition
%         if()
        f = f0;
%         end
% %            for tsteps = time
               f_eq = f_equilibrium_1d(rho_0,ux,v,t,RRR);
                % initialize variables
                 u_next = zeros(1,nx);
                 u_eq = zeros(1,nx);
                 u_l=zeros(1,nx);
                 u = zeros(1,nx);
                 for i = 1:nv
                      % load subcase
                      u_eq(:) = f_eq(i,:);
                      u(:) = f(i,:).*h;
                      u_l(:)=(1.-h).*f(i,:);
                      % Compute the smoothness factors, r(j), from data, u(j).
                       [r] = theta1d(u,a(i));
                        % Compute the Flux Limiter
                       [phi] = fluxlimiter1d(r,1); % using limiter = 1
                       % Compute TVD Fluxes
                       [F_left,F_right] = TVDflux1d(u,a(i),dtdx,phi);
                        [FL_left,FL_right] = TVDflux1d(u_l,a(i),dtdx,phi);
                       % Compute next time step
                        u_next= u - dtdx*(F_right - F_left).*h -dtdx*(FL_right - FL_left).*h + (dt/r_time).*h.*(u_eq-f(i,:));

                        % BC.4
                        u_next(1) = u_next(2);
                        u_next(nx) = u_next(nx-1);
                        % UPDATE info
                         u = u_next;

%                          for i=1:nx
%                              u
%                          end
                        % Going back to f
                         f(i,:) = u(:);
                 end
                   % Compute macroscopic moments

                   [n,j_x,E] = macromoments1d(k,w,f,v);

                   % UPDATE macroscopic properties
                     [ux,t,p,yun] = macroproperties1d(n,j_x,E,nx,nv,theta);
%                      [r,ux,t,p,yun] = macroproperties1d(n,j_x,E,nx,nv,theta);
%                         [p,yun] = macroproperties1d(n,j_x,E,nx,nv,theta);
%            end
%           case{2} % WENO k = 3 i.e. O(h^5)
%   otherwise
%          error('Order must be between 1 and 2');
% end


text:
function a = carry ( n, alpha )

%*****************************************************************************80
%
%% CARRY returns the CARRY matrix.
%
%  Discussion:
%
%    We assume that arithmetic is being done in base ALPHA.  We are adding
%    a column of N digits base ALPHA, as part of adding N random numbers.
%    We know the carry digit, between 0 and N-1, that is being carried into the
%    column sum (the incarry digit), and we want to know the probability of
%    the various carry digits 0 through N-1 (the outcarry digit) that could
%    be carried out of the column sum.
%
%    The carry matrix summarizes this data.  The entry A(I,J) represents
%    the probability that, given that the incarry digit is I-1, the
%    outcarry digit will be J-1.
%
%  Formula:
%
%    A(I,J) = ( 1 / ALPHA )^N * sum ( 0 <= K <= J-1 - floor ( I-1 / ALPHA ) )
%      (-1)^K * C(N+1,K) * C(N-I+(J-K)*ALPHA, N )
%
%  Example:
%
%    ALPHA = 10, N = 4
%
%    0.0715 0.5280 0.3795 0.0210
%    0.0495 0.4840 0.4335 0.0330
%    0.0330 0.4335 0.4840 0.0495
%    0.0210 0.3795 0.5280 0.0715
%
%  Square Properties:
%
%    A is generally not symmetric: A' /= A.
%
%    A is a Markov matrix.
%
%    A is centrosymmetric: A(I,J) = A(N+1-I,N+1-J).
%
%    LAMBDA(I) = 1 / ALPHA^(I-1)
%
%    det ( A ) = 1 / ALPHA^((N*(N-1))/2)
%
%    The eigenvectors do not depend on ALPHA.
%
%    A is generally not normal: A' * A /= A * A'.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 September 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    John Holte,
%    Carries, Combinatorics, and an Amazing Matrix,
%    The American Mathematical Monthly,
%    February 1997, pages 138-149.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer ALPHA, the numeric base being used in the addition.
%
%    Output, real A(N,N), the matrix.
%
  alpha = floor ( alpha );

  for i = 1 : n
    for j = 1 : n

      temp = 0.0;
      s = -1.0;

      for k = 0 : j - 1 - floor ( ( i - 1 ) / alpha )
        s = - s;
        c1 = r8_choose ( n + 1, k );
        c2 = r8_choose ( n - i + ( j - k ) * alpha, n );
        temp = temp + s * c1 * c2;
      end

      a(i,j) = temp / alpha^n;

    end

  end

  return
end


text:
function varargout = plot3(f, g, h, varargin)
%PLOT3   Plot for CLASSICFUN objects in 3-D space.
%   PLOT() is a three-dimensional analogue of PLOT().
%
%   PLOT3(X, Y, Z), where Z, Y, and Z are three CLASSICFUN objects, plots a line in
%   3-space. X, Y, and Z may be array-valued, but must have the same number of
%   columns.
%
%   Various line types, plot symbols, and colors may be obtained with PLOT3(X,
%   Y, Z, S) where s is a 1, 2 or 3 character string made from the characters
%   listed under the PLOT command.
%
%   H1 = PLOT3(X, Y, Z) returns a column vector of handles to lineseries
%   objects, one handle per plotted line (in the case of vector-valued CLASSICFUN
%   objects). [H1, H2] returns a second vector of column handles, this time for
%   each of the marker plots.
%
% See also PLOT, PLOTDATA.

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

% Deal with an empty input:
if ( isempty(f) )
    if ( nargout == 1 )
        varargout{1} = plot([]);
    end
    return
end

% Store the hold state of the current axis:
holdState = ishold;

%%
% We can only pl0t real data in this way:
f = real(f);
g = real(g);
h = real(h);

% Get the data for plotting from PLOTDATA():
data = plotData(f, g, h);

%%
% Plot the curve:
h1 = plot3(data.xLine, data.yLine, data.zLine, varargin{:});

% No points on the line plot:
set(h1, 'Marker', 'none')

% Hold the plot:
hold on

%%
% Plot the points:
h2 = plot3(data.xPoints, data.yPoints, data.zPoints, varargin{:});

% Change the style accordingly:
set(h2, 'LineStyle', 'none')
if ( all(strcmp(get(h2, 'Marker'), 'none')) && length(f) < 257 )
    set(h2,'Marker', 'o')
end

%%
% Return hold state to what it was before:
if ( ~holdState )
    hold off
end

% Give an output if one was requested:
if ( nargout > 0 )
    varargout{1} = h1;
    varargout{2} = h2;
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function item = cfg_entry(varargin)

% This is the entry configuration item class
%
% Data structure
% ==============
% Description fields
%    * name  - display name of config item
%    * tag   - tag of the menu item
%    * val   - 1x1 cell array
%    * check - (optional) function handle to implement configuration
%              specific subsasgn checks based on the harvested subtree
%              rooted at this node
%    * help  - help text
% GUI/job manager fields
%    * expanded
%    * hidden
% All fields above are inherited from the generic configuration item class.
%    * strtype
%    * num     - A 1-by-ndims vector of non-negative numbers, describing
%                the expected dimensions of the input. Dimensions with a
%                .num value of 0 or Inf can have an arbitrary number of
%                elements. In case of 2 dimensions, .val inputs will be
%                tried to match in un-transposed order first. If this
%                does not work, then .val inputs will be transposed and
%                matched again. If num is an empty matrix, no dimension
%                and size checks will be performed.
%                If strtype is 's' and num has 2 elements, these 2
%                elements code the min/max length of a string. This is a
%                workaround - in future versions num may be changed to a
%                2-by-ndims array encoding min/max values for each
%                dimension.
%    * def
%    * extras  - Extra information used for content checks of item.val{1}. The
%                following strtypes can use this extra information:
%                's' - a cell array of regular expressions. The val
%                string must match at least one of the regular
%                expressions.
%                'e' - a function handle [sts val] = f(val, item.num)
%
% Public Methods
% ==============
%    * get_strings - returns name of object
%    * gettag      - returns tag
%    * help        - returns help text
%    * harvest     - returns item.val{1}, or '<UNDEFINED>' if empty, see below
%    * all_set     - returns ~isempty(item.val), checks numel(item.val{1})
%                    against item.num
%
% Public internal Methods
% =======================
%    * subsasgn
%    * subsref
%
% 'strtype'
% The type of values that are entered by typing.  e.g. 'e' for evaluated.
% The valid value types are:
%   's'   string
%   'e'   evaluated expression - this can be any expression, even a struct
%          or cell
%   'f'   function or function handle
%   'n'   natural number (1..n)
%   'w'   whole number (0..n)
%   'i'   integer
%   'r'   real number
%   The following types are supported too, but there are no special
%   checks for validity of contents
%   'c'   indicator vector (e.g., 0101... or abab...)
%   'x'   contrast matrix
%   'p'   permutation
%
% Subscript Assignment Checks
% ===========================
% .num must conform to the semantics described above.
% The contents of item.val{1} will be checked to match .num and .strtype
% restrictions.
%
% Output in Job Structure (harvest)
% =================================
% cfg_entry uses cfg_item/harvest. If multiple dependencies are present
% and all can be resolved, the result will be a concatenation of all
% inputs. If concatenation fails, dependencies will not be resolved.
%
% The layout of the configuration tree and the types of configuration items
% have been kept compatible to a configuration system and job manager
% implementation in SPM5 (Statistical Parametric Mapping, Copyright (C)
% 2005 Wellcome Department of Imaging Neuroscience). This code has been
% completely rewritten based on an object oriented model of the
% configuration tree.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg_entry.m 5678 2013-10-11 14:58:04Z volkmar $

rev = '$Rev: 5678 $'; %#ok

myclass = mfilename;
% Get local fields and defaults from private/mysubs_fields
[fn, defs] = mysubs_fields;

if nargin == 1
    if isstruct(varargin{1})
        % assume input is a struct to be converted back into a class
        % return with error if this does not work
        if numel(fieldnames(varargin{1})) == numel(fn)+2 && ...
                all(isfield(varargin{1}, [fn(:)', {'cfg_item' 'cfg_leaf'}]))
            gitem = varargin{1}.cfg_item;
            sitem = rmfield(varargin{1},{'cfg_item', 'cfg_leaf'});
            item  = class(sitem, myclass, gitem, cfg_leaf);
            return;
        else
            cfg_message('matlabbatch:constructor:reclassify', ['Don''t know how to convert this ' ...
                            'into class ''%s''.'], myclass);
        end;
    end;
    if isa(varargin{1},myclass)
        item = varargin{1};
        return;
    end;
end;

mxpnargin = 4; % Max 4 arguments to parent initialisation
pnargin = min([nargin,mxpnargin]);
switch nargin
    case 0
        gitem = cfg_item;
    case {1,2,3,4,5,6,7,8,9}
        gitem = cfg_item(varargin{1:pnargin});
    otherwise
        cfg_message('matlabbatch:constructor:nargin', 'Wrong number of arguments.');
end;
for k=1:numel(fn)
    sitem.(fn{k})=defs{k};
end;
item = class(sitem, myclass, gitem, cfg_leaf);
if nargin > mxpnargin
    item.cfg_item.val{1} = varargin{mxpnargin+1};
    mxpnargin = mxpnargin+1;
end;
% set additional fields (if any) - field order as in mysubs_fields
for k = 1:min(numel(fn),nargin-mxpnargin)
    item.(fn{k}) = varargin{k+mxpnargin};
end;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% toy test
g1 = [ 1     2     3     1     2     3 ;
       4     5     6     5     6     4 ;
       3     3     3     1     1     5 ];
n = 6;
e1 = csaAssign(n,g1)

% big random test
n = 1000;
dg = (rand(n,n) > 0.5);
m = sum(dg(:));
i = find(dg==1)' - 1;
g2 = [ 1 + floor(i/n) ;
       1 + mod(i,n) + n ;
       1 + floor(rand(1,m)*1000) ];
tic;
e2 = csaAssign(2*n,g2);
toc;
if sum(e2(1,:)) ~= n*(n+1)/2, error('bug'); end
if sum(e2(2,:)) ~= n*(n+1)/2 + n*n, error('bug'); end
if sum(sum(e2(1:2,:))) ~= 2*n*(2*n+1)/2, error('bug'); end
disp('[n m cost] = ');
[n m sum(e2(3,:))]

author: Cloud-CV

repo: object-proposals


text:
function x = logistic_cdf_inv ( cdf, a, b )

%*****************************************************************************80
%
%% LOGISTIC_CDF_INV inverts the Logistic CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real CDF, the value of the CDF.
%    0.0 <= CDF <= 1.0.
%
%    Input, real A, B, the parameters of the PDF.
%    0.0 < B.
%
%    Output, real X, the corresponding argument.
%
  if ( cdf < 0.0 | 1.0 < cdf )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LOGISTIC_CDF_INV - Fatal error!\n' );
    fprintf ( 1, '  CDF < 0 or 1 < CDF.\n' );
    error ( 'LOGISTIC_CDF_INV - Fatal error!' );
  end

  x = a - b * log ( ( 1.0 - cdf ) / cdf );

  return
end


text:
function ex2bvp
%EX2BVP  Example 2 of the BVP tutorial.
%   A standard linear problem with a boundary layer at the origin.
%   The differential equation y'' + 3*p*y/(p + t^2)^2 = 0 has the
%   analytical solution y(t) = t/sqrt(p + t^2).  The parameter p
%   is taken to be 1e-5, a common value in tests.  The solution is
%   to have specified values at t = -0.1 and +0.1, values taken from
%   this analytical solution.
%
%   The default RelTol of 1e-3 gives an acceptable solution, but
%   reducing RelTol to 1e-4 resolves better the boundary layer.  A
%   constant guess is used for RelTol = 1e-3.  The same guess could be
%   used for RelTol = 1e-4, but a very much better guess is provided
%   by the solution previously computed for RelTol = 1e-3.

% Copyright 2002, The MathWorks, Inc.

% Evaluate the analytical solution for comparison.
tt = -0.1:0.01:+0.1;
p = 1e-5;
yy = tt ./ sqrt(p + tt .^2);

options = bvpset('stats','on','Fjacobian',@ex2Jac);

% BVPINT is used to specify an initial guess for the mesh of 10
% equally spaced points.  A constant guess based on a straight line
% between the boundary values for y is 0 for y(t) and 10 for y'(t).
solinit = bvpinit(linspace(-0.1,0.1,10),[0 10]);

sol = bvp4c(@ex2ode,@ex2bc,solinit, options);
t = sol.x;
y = sol.y;

figure
plot(t,y(1,:),tt,yy,'*')
axis([-0.1 0.1 -1.1 1.1])
title(['Linear boundary layer problem with RelTol = 1e-3.'])
xlabel('t')
ylabel('y and analytical (*) solutions')

fprintf('\n');

% A smaller RelTol is used to resolve better the boundary layer.
% The previous solution provides an excellent guess.
options = bvpset(options,'RelTol',1e-4);
sol = bvp4c(@ex2ode,@ex2bc,sol,options);
t = sol.x;
y = sol.y;

figure
plot(t,y(1,:),tt,yy,'*')
axis([-0.1 0.1 -1.1 1.1])
title(['Linear boundary layer problem with RelTol = 1e-4.'])
xlabel('t')
ylabel('y and analytical (*) solutions')

% --------------------------------------------------------------------------

function dydt = ex2ode(t,y)
%EX2ODE  ODE function for Example 2 of the BVP tutorial.
%   The components of y correspond to the original variables
%   as  y(1) = y, y(2) = y'.
p = 1e-5;
dydt = [ y(2)
        -3*p*y(1)/(p+t^2)^2];

% --------------------------------------------------------------------------

function dfdy = ex2Jac(t,y)
%EX2JAC  The Jacobian of the ODE function for Example 2 of the BVP tutorial.
p = 1e-5;
dfdy = [          0        1
         -3*p/(p+t^2)^2    0 ];

% --------------------------------------------------------------------------

function res = ex2bc(ya,yb)
%EX2BC  Boundary conditions for Example 2 of the BVP tutorial.
%   The boundary conditions are that the solution should agree
%   with the values of an analytical solution at both a and b.
p = 1e-5;
yatb = 0.1/sqrt(p + 0.01);
yata = - yatb;
res = [ ya(1) - yata
        yb(1) - yatb ];


text:
function lm = me_multin(Ns, epsil, A);
% Returns ldpc achevability data for a given (n,epsil,A)
cap = log2(1+A^2)/2;
Pdb = 20*log10(A);

lm = [];
for n = Ns;
	% these are absolute bounds on rate
	Rmin = 25/74; Rmax = 25/32;
	Rmax = min(Rmax, cap);
	% Step 1: find minimal acceptable rate
	Ptry = +Inf;
	Rnext = Rmin; Rstep = 0.1;
	while ~(Ptry <= Pdb)
		Rtry = Rnext;
		Ktry = floor(n*Rtry);
		[a b] = compute_ab(Ktry, Rtry);
		Ptry = (-norminv(epsil))/a + b;
		%% disp(sprintf('   Rtry = %.3g, Ktry = %d, Ptry = %.3g', Rtry, Ktry, Ptry));
		Rnext = Rnext + Rstep;
		if(Rnext > Rmax)
			if (Rstep == 0.1)
				Rnext = Rmin;
				Rstep = 0.01;
			else
				disp(sprintf(['me_multin(n = %d, epsil = %g, P = %.2g dB:' ...
						' could not find minrate'], n, epsil, Pdb));
				Rtry = Inf;
				break;
			end
		end
	end
	if(Rtry == Inf)
		lm = [lm NaN];
		continue;
	end
	Rmin = Rtry;
	Kmin = Ktry;
	Pmin = Ptry;
	% Step 2: find maximal acceptable rate

	Ptry = -Inf;
	Rnext = Rmax; Rstep = 0.1;
	while ~(Ptry >= Pdb)
		Rtry = Rnext;
		Ktry = ceil(n*Rtry);
		[a b] = compute_ab(Ktry, Rtry);
		Ptry = (-norminv(epsil))/a + b;
		Rnext = Rnext - Rstep;
		if(Rnext < Rmin)
			if (Rstep == 0.1)
				Rnext = Rmax;
				Rstep = 0.01;
			else
				disp(sprintf(['me_multin(n = %d, epsil = %g, P = %.2g dB:' ...
						' could not find maxrate'], n, epsil, Pdb));
				Rtry = Inf;
				break;
			end
		end
	end
	if(Rtry == Inf)
		lm = [lm NaN];
		continue;
	end
	Rmax = Rtry;
	Kmax = Ktry;
	Pmax = Ptry;
	% Step 3: proceed by binary division

	disp(sprintf(['me_multin(n = %d, eps = %g, P = %.2g dB: R = (%.3g, %.3g), '...
			'K = (%d, %d), P = (%.3g, %.3g)'], ...
			n, epsil, Pdb, Rmin, Rmax, Kmin, Kmax, Pmin, Pmax));

	while (Kmax - Kmin) > 1
		Rtry = (Rmin + Rmax)/2;
		Ktry = floor(Rtry*n);
		[a b] = compute_ab(Ktry, Rtry);
		Ptry = (-norminv(epsil))/a + b;
		%% disp(sprintf('   Rtry = %.3g, Ktry = %d, Ptry = %.3g', Rtry, Ktry, Ptry));
		if(Ptry == NaN) break; end
		if(Ptry > Pdb)
			Rmax = Rtry; Kmax = Ktry;
		else
			Rmin = Rtry; Kmin = Ktry;
		end
	end
	reason = 'PREC OK:';
	if(Ptry == NaN)
		reason = 'NAN BREAK:';
	end
	disp(sprintf('     %s lm = %d', reason, Kmin));
	lm = [lm Kmin];
end

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre


text:
classdef WOF_WeightIndividual < handle
% WOF_WeightIndividual - The class of an individual used in WOF to store
% weight variables.

% -----------------------------------------------------------------------
%  Copyright (C) 2020 Heiner Zille
%
%  This work is licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 4.0
%  International License. (CC BY-NC-SA 4.0). To view a copy of this license,
%  visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or see the
%  pdf-file "License-CC-BY-NC-SA-4.0.pdf" that came with this code.
%
%  You are free to:
%  * Share ? copy and redistribute the material in any medium or format
%  * Adapt ? remix, transform, and build upon the material
%  Under the following terms:
%  * Attribution ? You must give appropriate credit, provide a link to the
%     license, and indicate if changes were made. You may do so in any reasonable
%     manner, but not in any way that suggests the licensor endorses you or your use.
%  * NonCommercial ? You may not use the material for commercial purposes.
%  * ShareAlike ? If you remix, transform, or build upon the material, you must
%    distribute your contributions under the same license as the original.
%  * No additional restrictions ? You may not apply legal terms or technological
%    measures that legally restrict others from doing anything the license permits.
%
%  Author of this Code:
%   Heiner Zille <heiner.zille@ovgu.de> or <heiner.zille@gmail.com>
%
%  This code is based on the following publications:
%
%  1) Heiner Zille
%     "Large-scale Multi-objective Optimisation: New Approaches and a Classification of the State-of-the-Art"
%     PhD Thesis, Otto von Guericke University Magdeburg, 2019
%     http://dx.doi.org/10.25673/32063
%
%  2) Heiner Zille and Sanaz Mostaghim
%     "Comparison Study of Large-scale Optimisation Techniques on the LSMOP Benchmark Functions"
%     IEEE Symposium Series on Computational Intelligence (SSCI), IEEE, Honolulu, Hawaii, November 2017
%     https://ieeexplore.ieee.org/document/8280974
%
%  3) Heiner Zille, Hisao Ishibuchi, Sanaz Mostaghim and Yusuke Nojima
%     "A Framework for Large-scale Multi-objective Optimization based on Problem Transformation"
%     IEEE Transactions on Evolutionary Computation, Vol. 22, Issue 2, pp. 260-275, April 2018.
%     http://ieeexplore.ieee.org/document/7929324
%
%  4) Heiner Zille, Hisao Ishibuchi, Sanaz Mostaghim and Yusuke Nojima
%     "Weighted Optimization Framework for Large-scale Mullti-objective Optimization"
%     Genetic and Evolutionary Computation Conference (GECCO), ACM, Denver, USA, July 2016
%     http://dl.acm.org/citation.cfm?id=2908979
%
%  This file is intended to work with the PlatEMO framework version 2.5.
%  Date of publication of this code: 06.04.2020
%  Last Update of this code: 06.04.2020
%  A newer version of this algorithm may be available. Please contact the author
%  or see http://www.ci.ovgu.de/Research/Codes.html.
%
% The files may have been modified in Feb 2021 by the authors of the Platemo framework to work with the Platemo 3.0 release.
% -----------------------------------------------------------------------
% This file is derived from its original version containied in the PlatEMO
% framework.
% -----------------------------------------------------------------------

    properties(SetAccess = private)
        dec;        % Decision variables of the individual
        obj;        % Objective values of the individual
        con;        % Constraint values of the individual
        add;        % Additional properties of the individual
        ind;        % the actual individual to extract later
    end
    methods
        %% Constructor
        function obj = WOF_WeightIndividual(variables, GlobalDummy, addValues)

            if nargin > 0
                xPrimeVars = GlobalDummy.xPrime.dec;
                xPrimeSize = size(xPrimeVars,2);

                obj = WOF_WeightIndividual;

                % Set the infeasible decision variables to boundary values
                variables  = max(min(variables,GlobalDummy.upper),GlobalDummy.lower);

                x = WOF_transformationFunctionMatrixForm(xPrimeVars,variables(GlobalDummy.G),GlobalDummy.xPrimeupper,GlobalDummy.xPrimelower, GlobalDummy.psi);

                Problem = WOF_WeightIndividual.Current();
                obj.dec = variables;
                obj.ind = Problem.Evaluation(x);
                obj.obj = obj.ind.obj;
                obj.con = obj.ind.con;

                if nargin > 2
                    CallStack = dbstack();
                    Field     = CallStack(2).name;
                    obj.add.(Field) = addValues;
                end
            end

        end
        %% Get the matrix of decision variables of the population
        function value = decs(obj)
        %decs - Get the matrix of decision variables of the population
        %
        %   A = obj.decs returns the matrix of decision variables of the
        %   population obj, where obj is an array of INDIVIDUAL objects.

            value = cat(1,obj.dec);
        end
        %% Get the matrix of objective values of the population
        function value = objs(obj)
        %objs - Get the matrix of objective values of the population
        %
        %   A = obj.objs returns the matrix of objective values of the
        %   population obj, where obj is an array of INDIVIDUAL objects.

            value = cat(1,obj.obj);
        end
        %% Get the matrix of constraint values of the population
        function value = cons(obj)
        %cons - Get the matrix of constraint values of the population
        %
        %   A = obj.cons returns the matrix of constraint values of the
        %   population obj, where obj is an array of INDIVIDUAL objects.

            value = cat(1,obj.con);
        end
        %% Get the matrix of additional property of the population
        function value = adds(obj,addValue)
        %adds - Get the matrix of additional property values of the population
        %
        %   A = obj.adds(AddProper) returns the matrix of the values of the
        %   additional property of the INDIVIDUAL objects obj. The name of
        %   the additional property is same to the function name of the
        %   caller, that is, the values of one additional property of the
        %   individuals can only be obtained by the function which created
        %   them. If any individual in obj does not contain the specified
        %   additional property, assign it a default value specified in
        %   AddProper.

            CallStack = dbstack();
            Field     = CallStack(2).name;
            value     = zeros(length(obj),size(addValue,2));
            for i = 1 : length(obj)
                if ~isfield(obj(i).add,Field)
                    obj(i).add.(Field) = addValue(i,:);
                end
                value(i,:) = obj(i).add.(Field);
            end
        end
    end
    methods(Static, Sealed)
        function obj = Current(obj)
        %Current - Get or set the current PROBLEM object.

            persistent Problem;
            if nargin > 0
                Problem = obj;
            end
            if nargout > 0
                obj = Problem;
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
cd lib/mex_files;
make_mswindows
cd ../..
addpath ./demo
addpath ./lib/m_files
addpath ./lib/mex_files


text:
function res = qdist(a,b)
%QDIST        Implements  q(a,b)  metrical distance
%  name  qdist  is used to avoid ambiguities with variable  q
%
%     res = qdist(a,b)          max( abs(inf(a)-inf(b)) , abs(sup(a)-sup(b)) )
%
% for real intervals      max( abs(inf(a)-inf(b)) , abs(sup(a)-sup(b)) )
% for complex intervals   qdist(real(a),real(b)) + qdist(imag(a),imag(b))
%

% written  10/16/98     S.M. Rump
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  if ~isa(a,'intval')
    a = intval(a);
  end
  if ~isa(b,'intval')
    b = intval(b);
  end

  if ~a.complex & ~b.complex
    res = max( abs(a.inf-b.inf) , abs(a.sup-b.sup) ) ;
  else
    res = qdist(real(a),real(b)) + qdist(imag(a),imag(b)) ;
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function test08 ( sample_num )

%*****************************************************************************80
%
%% TEST08 times R4_EXP.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 May 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST08\n' );
  fprintf ( 1, '  Measure the time it takes R4_EXP to generate\n' );
  fprintf ( 1, '  %d exponential deviates.\n', sample_num );

  [ ke, fe, we ] = r4_exp_setup ( );

  seed = uint32 ( 123456789 );

  time1 = cputime;

  for sample = 1 : sample_num
    [ value, seed ] = r4_exp ( seed, ke, fe, we );
  end

  time2 = cputime;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  %f seconds.\n', time2 - time1 );

  return
end


text:
%  Figure 3.13      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
% script to generate Fig. 3.13
%  fig3_13.m
clf;
einv=1/exp(1);
num=1;
den=[1 1];
t=0:.05:4;
y=impulse(num,den,t);

% define some lines for the plot
tl=[0 1];
yl=[1 0];
t1=[1 1];
y1=[0 einv];
t2=[0 1];
y2=[einv einv];
figure();
plot(t,y,'-',tl,yl,'--',t1,y1,':',t2,y2,':','LineWidth',2)
title('Fig. 3.13(a) First order system impulse response')
xlabel('Time (sec)')
ylabel('h(t)')
text(0.7,0.6,'e^{-\sigmat}');
text(1.1,0.3679,'\leftarrow 1/e');
text(1,0.05,'\downarrow t= \tau');
% grid
nicegrid
pause;
% Figure 3.13 (b)
a=1;
num=[a];              % form numerator
den=[1 a];            % form denominator
t=0:0.01:4;           % form time vector
sys=tf(num,den);      % form system
h=impulse(sys,t);     % compute impulse response
figure();
plot(t,h);            % plot impulse response
y=step(sys,t);        % compute step response
hold
plot(t,y,'LineWidth',2);            % plot step response
xlabel('Time (sec)');
ylabel('h(t),y(t)');
title('Fig. 3.13(b) Impulse and step responses');
text(2,0.8,'y(t)');
text(2,0.2,'h(t)');
% grid
nicegrid