text:
% PART 1: Microchips classification ====================================

% Initialization
clear; close all; clc;

% Load the data --------------------------------------------------------
fprintf('Loading the data...\n\n');

data = load('microchips_tests.csv');

X = data(:, 1:2);
y = data(:, 3);

% Plotting the data ----------------------------------------------------
fprintf('Plotting the data...\n\n');

% Find indices of positive and negative examples.
positiveIndices = find(y == 1);
negativeIndices = find(y == 0);

% Plot examples.
hold on;
plot(X(positiveIndices, 1), X(positiveIndices, 2), 'k+', 'LineWidth', 2, 'MarkerSize', 7);
plot(X(negativeIndices, 1), X(negativeIndices, 2), 'ko', 'MarkerFaceColor', 'y', 'MarkerSize', 7);

% Draw labels and Legend
xlabel('Microchip Test 1');
ylabel('Microchip Test 2');
legend('y = 1', 'y = 0');

% Running logistic regression ------------------------------------------
fprintf('Running logistic regression...\n\n');

% Add more polynomial features in order to make decision boundary to have more complex curve form.
polynomial_degree = 6;
X = add_polynomial_features(X(:, 1), X(:, 2), polynomial_degree);

% Run the regression.
lambda = 1;
[theta, J, J_history, exit_flag] = logistic_regression_train(X, y, lambda);

fprintf('- Initial cost: %f\n', J_history(1));
fprintf('- Optimized cost: %f\n\n', J);

% Plotting decision boundaries -----------------------------------------
fprintf('Plotting decision boundaries...\n\n');

% Generate a grid range.
u = linspace(-1, 1, 50);
v = linspace(-1, 1, 50);
z = zeros(length(u), length(v));
% Evaluate z = (x * theta) over the grid.
for i = 1:length(u)
    for j = 1:length(v)
        % Add polinomials.
        x = add_polynomial_features(u(i), v(j), polynomial_degree);
        % Add ones.
        x = [ones(size(x, 1), 1), x];
        z(i, j) = x * theta;
    end
end

% It is mportant to transpose z before calling the contour.
z = z';

% Plot z = 0
% Notice you need to specify the range [0, 0]
contour(u, v, z, [0, 0], 'LineWidth', 2);
title(sprintf('lambda = %g', lambda));
legend('y = 1', 'y = 0', 'Decision boundary');

hold off;

% Trying to predict custom experiments ------------------------------------
fprintf('Trying to predict custom experiments...\n\n');

x = [
    0, 0;
    -0.5, -0.5
];

% Add polinomials.
x = add_polynomial_features(x(:, 1), x(:, 2), polynomial_degree);
% Add ones.
x = [ones(size(x, 1), 1), x];

probabilities = hypothesis(x, theta);
fprintf(' %f \n', probabilities);

fprintf('\nProgram paused. Press enter to train logistic regression to recognize digits.\n');
pause;

% PART 2: Handwritten digits classification ==============================
clear; close all; clc;

% Load training data
fprintf('Loading training data...\n');
load('digits.mat');

% Plotting some training example ----------------------------------------------------
fprintf('Visualizing data...\n');

% Randomly select 100 data points to display
random_digits_indices = randperm(size(X, 1));
random_digits_indices = random_digits_indices(1:100);

display_data(X(random_digits_indices, :));

% Setup the parameters you will use for this part of the exercise
input_layer_size = 400;  % 20x20 input images of digits.
num_labels = 10; % 10 labels, from 1 to 10 (note that we have mapped "0" to label 10).

fprintf('Training One-vs-All Logistic Regression...\n')
lambda = 0.01;
num_iterations = 50;
[all_theta] = one_vs_all(X, y, num_labels, lambda, num_iterations);

fprintf('Predict for One-Vs-All...\n')
pred = one_vs_all_predict(all_theta, X);

fprintf('\nTraining Set Accuracy: %f\n', mean(double(pred == y)) * 100);

author: trekhleb

repo: machine-learning-octave


text:
classdef CEC2020_F3 < PROBLEM
% <single> <real>
% Shifted and rotated Lunacek bi-Rastrigin function

%------------------------------- Reference --------------------------------
% C .T. Yue, K. V. Price, P. N. Suganthan, J. J. Liang, M. Z. Ali, B. Y.
% Qu, N. H. Awad, and P. P Biswas, Problem definitions and evaluation
% criteria for the CEC 2020 special session and competition on single
% objective bound constrained numerical optimization, Zhengzhou University,
% China and Nanyang Technological University, Singapore, 2019.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties
        O;      % Optimal decision vector
        Mat;	% Rotation matrix
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            CallStack = dbstack('-completenames');
            load(fullfile(fileparts(CallStack(1).file),'CEC2020.mat'),'Data');
            obj.O = Data{3}.o;
            obj.M = 1;
            if isempty(obj.D) || obj.D < 10
                obj.D   = 5;
                obj.Mat = Data{3}.M_5;
            elseif obj.D < 15
                obj.D   = 10;
                obj.Mat = Data{3}.M_10;
            elseif obj.D < 20
                obj.D   = 15;
                obj.Mat = Data{3}.M_15;
            else
                obj.D   = 20;
                obj.Mat = Data{3}.M_20;
            end
            obj.lower    = zeros(1,obj.D) - 100;
            obj.upper    = zeros(1,obj.D) + 100;
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            s   = 1 - 1/(2*sqrt(obj.D+20)-8.2);
            mu0 = 2.5;
            mu1 = -sqrt((mu0^2-1)/s);
            Y   = (PopDec-repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1))/10;
            tmp = 2*repmat(sign(obj.O(1:size(PopDec,2))),size(PopDec,1),1).*Y + mu0;
            Z   = (tmp-mu0)*obj.Mat';
            PopObj = 700 + min(sum((tmp-mu0).^2,2),obj.D+s*sum((tmp-mu1).^2,2)) + 10*(obj.D-sum(cos(2*pi*Z),2));
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function varargout = gmm(varargin)
% Vl_GMM  Learn a Gaussian Mixture Model using EM
%   [MEANS, COVARIANCES, PRIORS] = VL_GMM(X, NUMCLUSTERS) fits a GMM with
%   NUMCLUSTERS components to the data X. Each column of X represent a
%   sample point. X may be either SINGLE or DOUBLE. MEANS, COVARIANCES, and
%   PRIORS are respectively the means, the diagonal covariances, and
%   the prior probabilities of the Guassian modes. MEANS and COVARIANCES
%   have the same number of rows as X and NUMCLUSTERS columns with one
%   column per mode. PRIORS is a row vector with NUMCLUSTER entries
%   summing to one.
%
%   [MEANS, COVARIANCES, PRIORS, LL] = VL_GMM(...) returns the
%   loglikelihood (LL) of the model as well.
%
%   [MEANS, COVARIANCES, PRIORS, LL, POSTERIORS] = VL_GMM(...) returns
%   the posterior probabilities POSTERIORS of the Gaussian modes given
%   each data point. The POSTERIORS matrix has NUMCLUSTERS rows and
%   NUMDATA columns.
%
%   VL_GMM() supports different initialization and optimization
%   methods. Specifically, the following options are supported:
%
%   Verbose::
%     Increase the verbosity level (may be specified multiple times).
%
%   Initialization:: RAND
%     RAND initializes the means as random data poitns and the
%     covaraince matrices as the covariance of X. CUSTOM allow
%     specifying the initial means, covariances, and prior
%     probabilities.
%
%   InitMeans:: none
%     Specify the initial means (size(X,1)-by-NUMCLUSTERS matrix).
%
%   InitPriors:: none
%     Specify the initial weights (a vector of dimension NUMCLUSTER).
%
%   InitCovariances:: none
%     Specify the initial diagonal covariance matrices
%
%   NumRepetitions:: 1
%     Number of times to restart EM. The solution with maximum
%     loglikelihood is returned.
%
%   CovarianceBound:: 10e-6
%     Set the lower bound on the diagonal covariance values.
%     The bound can be either a scalar or a vector with one
%     entry per dimension. Using null bounds is possible, but
%     may yield degenerate solutions, including NaNs.
%
%   Example::
%     VL_GMM(X, 10, 'verbose', 'MaxNumIterations', 20) learns a
%     mixture of ten Gaussians using at most twenty iterations of the
%     algorithm.
%
%   See also: GMMs, VL_KMEANS(), VL_HELP().
[varargout{1:nargout}] = vl_gmm(varargin{:});

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function [rhsu] = HeatLDGRHS1D(u,time)

% function [rhsu] = HeatLDGRHS1D(u,time,a,ax)
% Purpose  : Evaluate RHS flux in 1D heat equation using an LDG flux

Globals1D;

% Define field differences at faces
du = zeros(Nfp*Nfaces,K); du(:) = (1.0+nx(:)).*(u(vmapM)-u(vmapP))/2.0;

% impose boundary condition -- Dirichlet BC's
%uin  = -u(vmapI); du(mapI) = (1.0+nx(mapI)).*(u(vmapI)- uin)/2.0;
%uout = -u(vmapO); du(mapO) = (1.0+nx(mapO)).*(u(vmapO)-uout)/2.0;
uin  = 0.1; du(mapI) = (1.0+nx(mapI)).*(u(vmapI)- uin)/2.0;
uout = 0.1; du(mapO) = (1.0+nx(mapO)).*(u(vmapO)-uout)/2.0;

% Compute q
q = rx.*(Dr*u) - LIFT*(Fscale.*(nx.*du));
dq = zeros(Nfp*Nfaces,K); dq(:) = (1.0-nx(:)).*(q(vmapM)-q(vmapP))/2.0;

% impose boundary condition -- Neumann BC's
%qin  = q(vmapI); dq(mapI) = (1.0-nx(mapI)).*(q(vmapI)- qin)/2.0;
%qout = q(vmapO); dq(mapO) = (1.0-nx(mapO)).*(q(vmapO)-qout)/2.0;
dq(mapI)=0;
dq(mapO)=0;

% compute right hand sides of the semi-discrete PDE
rhsu = rx.*(Dr*q) - LIFT*(Fscale.*(nx.*dq));
return

repo: Aero-matlab


text:
function plot_llkstest(a,time,timef,bootloops,maepi)
    % function plot_llkstes(a,time,timef,bootloops,maepi);
    % --------------------------------------------------
    % Plots Ncum observed vs. Ncum modeled for specified time windows
    % with choosing the model for the learning period and performs a KS-Test
    %
    % Input variables:
    % a         : earthquake catalog
    % time      : learning period fo fit Omori parameters
    % timef     : forecast period
    % bootloops : Number of bootstraps
    % maepi     : mainshock
    %
    % J.Woessner
    % last update: 20.07.04

    % Surpress warnings from fmincon
    warning off;

    %[m_main, main] = max(a.Magnitude);
    date_matlab = datenum(a.Date.Year,a.Date.Month,a.Date.Day,a.Date.Hour,a.Date.Minute,zeros(size(a,1),1));
    date_main = datenum(floor(maepi(3)),maepi(4),maepi(5),maepi(8),maepi(9),0);
    time_aftershock = date_matlab-date_main;

    % Aftershock catalog
    vSel1 = time_aftershock(:) > 0;
    tas = time_aftershock(vSel1);
    eqcatalogue = a.subset(vSel1);

    % Estimation of Omori parameters from learning period
    l = tas <= time;
    time_as=tas(l);
    % Times up to the forecast time
    lf = tas <= time+timef ;
    time_asf= [tas(lf) ];
    time_asf=sort(time_asf);

    % Select biggest aftershock earliest in time, but more than 1 day after
    % mainshock and in learning period
    mAfLearnCat = eqcatalogue(l,:);
    fDay = 1;
    ft_c=fDay/365; % Time not considered to find biggest aftershock
    vSel = (mAfLearnCat(:,3) > maepi(:,3)+ft_c & mAfLearnCat(:,3)<= maepi(:,3)+time/365);
    mCat = mAfLearnCat(vSel,:);
    vSel = mCat(:,6) == max(mCat(:,6));
    vBigAf = mCat(vSel,:);
    if length(mCat(:,1)) > 1
        vSel = vBigAf(:,3) == min(vBigAf(:,3));
        vBigAf = vBigAf(vSel,:);
    end

    date_biga = datenum(floor(vBigAf(3)),vBigAf(4),vBigAf(5),vBigAf(8),vBigAf(9),0);
    fT1 = date_biga - date_main; % Time of big aftershock

    % Calculate p,c,k for dataset
    prompt  = {'Enter model number (1:pck, 2:pckk, 3:ppckk, 4:ppcckk:'};
    title   = 'Model selection for fitting aftershock sequence';
    lines= 1;
    def     = {'1'};
    answer  = inputdlg(prompt,title,lines,def);
    nMod = str2double(answer{1});

    % Calculate uncertainty and mean values of p,c,and k
    [mMedModF, mStdL, loopout] = brutebootloglike_a2(time_as, time_asf, bootloops,fT1,nMod);
    pmed1 = mMedModF(1,1);
    pmed2 = mMedModF(1,3);
    cmed1 = mMedModF(1,5);
    cmed2 = mMedModF(1,7);
    kmed1 = mMedModF(1,9);
    kmed2 = mMedModF(1,11);

    % Compute model according to model choice
    if nMod == 1
        [pval1, pval2, cval1, cval2, kval1, kval2, fAIC, fL] = bruteforceloglike_a2(time_as,fT1,nMod);
    elseif nMod == 2
        [pval1, pval2, cval1, cval2, kval1, kval2, fAIC, fL] = bruteforceloglike_a2(time_as,fT1,nMod);
    elseif nMod == 3
        [pval1, pval2, cval1, cval2, kval1, kval2 , fAIC, fL] = bruteforceloglike_a2(time_as,fT1,nMod);
    else
        [pval1, pval2, cval1, cval2, kval1, kval2, fAIC, fL] = bruteforceloglike_a2(time_as,fT1,nMod);
    end

    % Start plotting
    figure_w_normalized_uicontrolunits('Numbertitle','off','Name','Aftershock modelling fit')

    % Plot the forecast ...
    cumnrf = (1:length(time_as))';
    cumnr_modelf = [];
    if nMod == 1
        for i=1:length(time_as)
            if pval1 ~= 1
                cm = kval1/(pval1-1)*(cval1^(1-pval1)-(time_as(i)+cval1)^(1-pval1));
            else
                cm = kval1*log(time_as(i)/cval1+1);
            end
            cumnr_modelf = [cumnr_modelf; cm];
        end % END of FOR on length(time_as)
    else
        for i=1:length(time_as)
            if time_as(i) <= fT1
                if pval1 ~= 1
                    cm = kval1/(pval1-1)*(cval1^(1-pval1)-(time_as(i)+cval1)^(1-pval1));
                else
                    cm = kval1*log(time_as(i)/cval1+1);
                end
                cumnr_modelf = [cumnr_modelf; cm];
            else
                if (pval1 ~= 1 & pval2 ~= 1)
                    cm = kval1/(pval1-1)*(cval1^(1-pval1)-(time_as(i)+cval1)^(1-pval1))+ kval2/(pval2-1)*(cval2^(1-pval2)-(time_as(i)-fT1+cval2)^(1-pval2));
                else
                    cm = kval1*log(time_as(i)/cval1+1) + kval2*log((time_as(i)-fT1)/cval2+1);
                end
                cumnr_modelf = [cumnr_modelf; cm];
            end %END of IF on fT1
        end % End of FOR length(time_as)
    end % End of if on nMod
    time_as=sort(time_as);
    cumnr_modelf=sort(cumnr_modelf);

    pf1 =  plot(time_as,cumnr_modelf,'g-.','Linewidth',2);
    hold on
    pf2 =  plot(time_as,cumnrf, 'b-','Linewidth',2);
    paf = plot(fT1, 0,'h','MarkerFaceColor',[1 1 0],'MarkerSize',12,'MarkerEdgeColor',[0 0 0] );

    % Calculate KSTEST2 as a measure of the goodness of fit
    [H,P,KSSTAT] = kstest2(cumnr_modelf,cumnrf)

    % Calculate RMS
    i=(1:1:length(time_as))';
    fRMS = (sum((i-cumnr_modelf).^2)/length(i))^0.5;

    % Round values for output
    pval1 = round(100*pval1)/100;
    pval2 = round(100*pval2)/100;
    cval1 = round(1000*cval1)/1000;
    cval2 = round(1000*cval2)/1000;
    kval1 = round(10*kval1)/10;
    kval2 = round(10*kval2)/10;
    pmed1 = round(100*pmed1)/100; mStdL(1,1) = round(100*mStdL(1,1))/100;
    pmed2 = round(100*pmed2)/100; mStdL(1,2) = round(100*mStdL(1,2))/100;
    cmed1 = round(1000*cmed1)/1000; mStdL(1,3) = round(1000*mStdL(1,3))/1000;
    cmed2 = round(1000*cmed2)/1000; mStdL(1,4) = round(1000*mStdL(1,4))/1000;
    kmed1 = round(10*kmed1)/10; mStdL(1,5) = round(100*mStdL(1,5))/100;
    kmed2 = round(10*kmed2)/10; mStdL(1,6)= round(100*mStdL(1,6))/100;
    fRMS = round(100*fRMS)/100;

    % Get Y limits for positioning texts
    yy = get(gca,'ylim');

    if nMod == 1
        string1=['p = ' num2str(pval1) '; c = ' num2str(cval1) '; k = ' num2str(kval1) ];
        string3=['pm = ' num2str(pmed1) '+-' num2str(mStdL(1,1)) '; cm = ' num2str(cmed1) '+-' num2str(mStdL(1,3)) '; km = ' num2str(kmed1) '+-' num2str(mStdL(1,5))];
        text(max(time_asf)*0.05,yy(2)*0.9,string1,'FontSize',10);
        text(max(time_asf)*0.05,yy(2)*0.8,string3,'FontSize',10);
    elseif nMod == 2
        string1=['p = ' num2str(pval1) '; c = ' num2str(cval1) '; k1 = ' num2str(kval1) '; k2 = ' num2str(kval2) ];
        string3=['pm = ' num2str(pmed1) '+-' num2str(mStdL(1,1)) '; cm = ' num2str(cmed1) '+-' num2str(mStdL(1,3)) '; km1 = ' num2str(kmed1) '+-' num2str(mStdL(1,5)) '; km2 = ' num2str(kmed2) '+-' num2str(mStdL(1,6))];
        text(max(time_asf)*0.05,yy(2)*0.9,string1,'FontSize',10);
        text(max(time_asf)*0.05,yy(2)*0.8,string3,'FontSize',10);
    elseif nMod == 3
        string1=['p1 = ' num2str(pval1) '; c = ' num2str(cval1) '; k1 = ' num2str(kval1) ];
        string2=['p2 = ' num2str(pval2) '; k2 = ' num2str(kval2) ];
        string3=['pm1 = ' num2str(pmed1) '+-' num2str(mStdL(1,1)) '; cm = ' num2str(cmed1) '+-' num2str(mStdL(1,3)) '; km1 = ' num2str(kmed1) '+-' num2str(mStdL(1,5))];
        string4=['pm2 = ' num2str(pmed2) '+-' num2str(mStdL(1,2)) '; km2 = ' num2str(kmed2) '+-' num2str(mStdL(1,6))];
        text(max(time_asf)*0.05,yy(2)*0.9,string1,'FontSize',10);
        text(max(time_asf)*0.05,yy(2)*0.85,string2,'FontSize',10);
        text(max(time_asf)*0.05,yy(2)*0.8,string3,'FontSize',10);
        text(max(time_asf)*0.05,yy(2)*0.75,string4,'FontSize',10);
    else
        string1=['p1 = ' num2str(pval1) '; c1 = ' num2str(cval1) '; k1 = ' num2str(kval1) ];
        string2=['p2 = ' num2str(pval2) '; c2 = ' num2str(cval2) '; k2 = ' num2str(kval2) ];
        string3=['pm1 = ' num2str(pmed1) '+-' num2str(mStdL(1,1)) '; cm1 = ' num2str(cmed1) '+-' num2str(mStdL(1,3)) '; km1 = ' num2str(kmed1) '+-' num2str(mStdL(1,5))];
        string4=['pm2 = ' num2str(pmed2) '+-' num2str(mStdL(1,2)) '; cm2 = ' num2str(cmed2) '+-' num2str(mStdL(1,4)) '; km2 = ' num2str(kmed2) '+-' num2str(mStdL(1,6))];
        text(max(time_asf)*0.05,yy(2)*0.9,string1,'FontSize',10);
        text(max(time_asf)*0.05,yy(2)*0.85,string2,'FontSize',10);
        text(max(time_asf)*0.05,yy(2)*0.8,string3,'FontSize',10);
        text(max(time_asf)*0.05,yy(2)*0.75,string4,'FontSize',10);
    end
    string=['H = ' num2str(H) ' P = ' num2str(P) ' KS-Statistic) = ' num2str(KSSTAT)];
    text(max(time_asf)*0.05,yy(2)*0.1,string,'FontSize',10);
    sAIC = ['AIC = ' num2str(fAIC)];
    text(max(time_asf)*0.05,yy(2)*0.05,sAIC,'FontSize',10);
    sRMS = ['RMS = ' num2str(fRMS)];
    text(max(time_asf)*0.05,yy(2)*0.15,sRMS,'FontSize',10);
    % Legend
    sModel = ['Model ' num2str(nMod)];
    legend([pf2 pf1 paf],'Data',sModel,'Sec. AF','location','Best');

    % Calculate fits of different models
    mRes = [];
    % Modified Omori law (pck)
    nMod = 1; [pval1, pval2, cval1, cval2, kval1, kval2, fAIC, fL] = bruteforceloglike_a2(time_as,fT1,nMod);
    mRes = [mRes; nMod, pval1, pval2, cval1, cval2, kval1, kval2, fAIC, fL];
    % MOL with secondary aftershock (pckk)
    nMod = 2; [pval1, pval2, cval1, cval2, kval1, kval2, fAIC, fL] = bruteforceloglike_a2(time_as,fT1,nMod);
    mRes = [mRes; nMod, pval1, pval2, cval1, cval2, kval1, kval2, fAIC, fL];
    % MOL with secondary aftershock (ppckk)
    nMod = 3; [pval1, pval2, cval1, cval2, kval1, kval2, fAIC, fL] = bruteforceloglike_a2(time_as,fT1,nMod);
    mRes = [mRes; nMod, pval1, pval2, cval1, cval2, kval1, kval2, fAIC, fL];
    % MOL with secondary aftershock (ppcckk)
    nMod = 4; [pval1, pval2, cval1, cval2, kval1, kval2, fAIC, fL] = bruteforceloglike_a2(time_as,fT1,nMod);
    mRes = [mRes; nMod, pval1, pval2, cval1, cval2, kval1, kval2, fAIC, fL];

    % Select best fitting model by AIC
    vSel = (mRes(:,8)==min(mRes(:,8)));
    mRes = mRes(vSel,:);
    if length(mRes(:,1)) > 1
        vSel = (mRes(:,1)==min(mRes(:,1)));
        mRes = mRes(vSel,:);
    end
    % Model to use for bootstrapping as of lowest AIC to observed data
    nMod = mRes(1,1);

    sModel1 = ['Info: Best model is ' num2str(nMod)];
    text(max(time_asf)*0.05,yy(2)*0.2,sModel1,'FontSize',10);
    % Figure settings
    set(gca,'FontSize',12,'Fontweight','bold','Linewidth',2)
    xlabel('Time [Days after mainshock]','FontSize',12,'Fontweight','bold')

author: CelsoReyes


text:
function eigg_segmentation
% This code was written by Muhammet Balcilar , France,
% muhammetbalcilar@gmail.com, inspired from following reference
%
%  This Matlab code implements an edge-based active contour model as an
%  application of the Distance Regularized Level Set Evolution (DRLSE) formulation in Li et al's paper:
%
%      C. Li, C. Xu, C. Gui, M. D. Fox, "Distance Regularized Level Set Evolution and Its Application to Image Segmentation",
%        IEEE Trans. Image Processing, vol. 19 (12), pp.3243-3254, 2010.

close all

%% step1, read grayscale image
Img=imread('Inputs/eigg_scotland.jpg');
Img=rgb2gray(Img);
frm=0;
%% step2, set params
timestep=5;  % time step
mu=0.2;  % coefficient of the distance regularization term R(phi)
lambda=5; %coefficient of the weighted length term L(phi)
alfa= -3;  %  coefficient of the weighted area term A(phi)
epsilon=1.5; % papramater that specifies the width of the DiracDelta function
c0=2;
maxiter=602;
sigma=2.0;    % scale parameter in Gaussian kernel

%% step3 smooth image with gaussian filter
G=fspecial('gaussian',30,sigma); % 15 Caussian kernel
Img_smooth=conv2(Img,G,'same');  % smooth image by Gaussiin convolution
figure(1);
imagesc(Img_smooth,[0, 255]); axis off; axis equal; colormap(gray);
title('Smoothed image');

%% step4 calculate edge indicator according to Eq23
[Ix,Iy]=gradient(Img_smooth);
f=Ix.^2+Iy.^2;
g=1./(1+f);  % edge indicator function.
g=exp(-f);
figure(2);
imagesc(g); axis off; axis equal;
title('g,  edge indicator');

%% step5, set initial phi
phi = -c0*ones(size(Img));
phi(50:650,50:500)=c0;
%phi(350:370,200:220)=-c0;
figure(3);
imagesc(phi);
axis off; axis equal;colormap(jet);
title('initial phi matrix');

[vx, vy]=gradient(g);
figure(4);
subplot(1,2,1);imagesc(vx); title('x directioned gradient of g');
subplot(1,2,2);imagesc(vy); title('y directioned gradient of g');

for k=1:maxiter
    %% step6, check boundary conditions
    phi=NeumannBoundCond(phi);

    %% step 7 calculate differential of regularized term in Eq.30
    distRegTerm=distReg_p2(phi);

    %% step8 calculate differential of area term in Eq.30
    diracPhi=Dirac(phi,epsilon);
    areaTerm=diracPhi.*g;

    %% step9 calculate differential of length term in Eq.30
    [phi_x,phi_y]=gradient(phi);
    s=sqrt(phi_x.^2 + phi_y.^2);
    Nx=phi_x./(s+1e-10); % add a small positive number to avoid division by zero
    Ny=phi_y./(s+1e-10);
    edgeTerm=diracPhi.*(vx.*Nx+vy.*Ny) + diracPhi.*g.*div(Nx,Ny);

    %% step 10 update phi according to Eq.20
    phi=phi + timestep*(mu/timestep*distRegTerm + lambda*edgeTerm + alfa*areaTerm);

     %% show result in every 50 iteration
    if mod(k,50)==1
        frm=frm+1;
        close all
        h=figure(5);
        set(gcf,'color','w');
        subplot(1,2,1);
        II=Img;
        II(:,:,2)=Img;II(:,:,3)=Img;
        imshow(II); axis off; axis equal; hold on;
        q=contour(phi, [0,0], 'r');
        msg=['contour result , iteration number=' num2str(k)];
        title(msg);
        subplot(1,2,2);
        mesh(-phi,[-2 2]);
        hold on;  contour(phi, [0,0], 'r','LineWidth',2);

        view([180-30 -65-180]);
        msg=['phi result , iteration number=' num2str(k)];
        title(msg);
        pause(0.1)

        % gif video

        frame = getframe(h);
        im = frame2im(frame);
        [imind,cm] = rgb2ind(im,256);
        %Write to the GIF File
        if frm == 1
            imwrite(imind,cm,'Outputs/eigg.gif','gif', 'Loopcount',inf);
        else
            imwrite(imind,cm,'Outputs/eigg.gif','gif','WriteMode','append');
        end
    end

    %% step 11 if maxiter done then finish, else return step6
end
%% Step 12. show last iteration results
figure(6);
imagesc(Img,[0, 255]); axis off; axis equal; colormap(gray); hold on;  contour(phi, [0,0], 'r');
msg=['phi result , iteration number=' num2str(k)];
title(msg);

function f = distReg_p2(phi)
% compute the distance regularization term with the double-well potential p2 in eqaution (16)
[phi_x,phi_y]=gradient(phi);
s=sqrt(phi_x.^2 + phi_y.^2);
a=(s>=0) & (s<=1);
b=(s>1);
ps=a.*sin(2*pi*s)/(2*pi)+b.*(s-1);  % compute first order derivative of the double-well potential p2 in eqaution (16)
dps=((ps~=0).*ps+(ps==0))./((s~=0).*s+(s==0));  % compute d_p(s)=p'(s)/s in equation (10). As s-->0, we have d_p(s)-->1 according to equation (18)
f = div(dps.*phi_x - phi_x, dps.*phi_y - phi_y) + 4*del2(phi);

function f = div(nx,ny)
[nxx,junk]=gradient(nx);
[junk,nyy]=gradient(ny);
f=nxx+nyy;

function f = Dirac(x, sigma)
f=(1/2/sigma)*(1+cos(pi*x/sigma));
b = (x<=sigma) & (x>=-sigma);
f = f.*b;

function g = NeumannBoundCond(f)
% Make a function satisfy Neumann boundary condition
[nrow,ncol] = size(f);
g = f;
g([1 nrow],[1 ncol]) = g([3 nrow-2],[3 ncol-2]);
g([1 nrow],2:end-1) = g([3 nrow-2],2:end-1);
g(2:end-1,[1 ncol]) = g(2:end-1,[3 ncol-2]);

author: balcilar

repo: DRLSE-Image-Segmentation


text:
classdef PartialDCT < spx.dict.Operator

properties(SetAccess=private)
    % Dimensions
    row_pics
    col_perm
    num_rows
    num_cols
end

methods

function self = PartialDCT(row_pics, col_perm)
    if nargin < 2
        error('Row selection and column permutation must be specified');
    end
    self.row_pics = row_pics;
    self.col_perm = col_perm;
    self.num_rows = numel(row_pics);
    self.num_cols = numel(col_perm);
    if self.num_rows > self.num_cols
        error('Number of rows cannot be greater than number of columns.');
    end
end

function [mm, nn]  = get_size(self)
    mm = self.num_rows;
    nn = self.num_cols;
end

function result = apply(self, vectors)
    if size(vectors, 1) ~= self.num_cols
        error('Dimensions mismatch');
    end
    result = dct(vectors(self.col_perm, :));
    result = result(self.row_pics, :);
end

function result = apply_ctranspose(self, vectors)
    if size(vectors, 1) ~= self.num_rows
        error('Dimensions mismatch');
    end
    num_vectors = size(vectors, 2);
    full_vectors = zeros(self.num_cols, num_vectors);
    full_vectors(self.row_pics, :) = vectors;
    result = zeros(self.num_cols, num_vectors);
    result(self.col_perm, :) = idct(full_vectors);
end

function result = double(self)
    % Converts the operator into a MATRIX
    result = self.apply(eye(self.num_cols));
end

function result = norm(self)
    error('Not implemented');
end

end

end

author: indigits

repo: sparse-plex


text:
%CellCycle medium test
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 11/23/2010
classdef CellCycle_Test < TestCase
    properties
        simulation
    end

    methods
        function this = CellCycle_Test(name)
            this = this@TestCase(name);
        end

        function setUp(this)
            this.simulation = edu.stanford.covert.cell.sim.SimulationFixture.load();
        end

        %Test that the expected lengths of the 3 phases (replication initation,
        %replication, cytokinesis) of the cell cycle sum to the expected cell
        %cycle length.
        %
        %The expected lengths of the 3 phases of the cell cycle have been
        %constructed as follows:
        %1. The expected length of the replication phase was calculated
        %   analytically using experimentally observed values for the
        %   replication constants.
        %2. The expected length of the cytokinesis phase was estimated by
        %   simulation. This is implemented in the cytokinesis medium test.
        %3. The expected length of the replication initiation phase was set to
        %   be the difference of the total cell cycle length and that of the
        %   replication and cytokinesis phases, and the replication initation
        %   constants (in particular the site and state cooperativity constants)
        %   were fit to achieve this expected replication initation length. This
        %   is implemented in the replication initation process test class.
        function testCellCycleLength(this)
            sim = this.simulation;

            %% lengths of individual phases of cell cycle
            %1. Replication initiation: replication initiation test asserts that
            %   replication initiation takes the expected time
            %2. Replication: replication unit test asserts that replication
            %   takes the expected time
            %3. Cytokinesis: Cytokinesis medium test asserts that cytokinesis
            %   takes the expected time

            %% total length of cell cycle
            time = sim.state('Time');
            assertEqual(time.cellCycleLength, ...
                + time.replicationInitiationDuration + ...
                + time.replicationDuration + ...
                + time.cytokinesisDuration);
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function  mShed = BSLshedMask(Wsheds,Cut)
%"BSLshedMask"
%   BSL subfunction -- Returns WS masks used for a BSL estimate made with a
%   user specified number of low activity regions removed
%
% CRS 05/20/13
%
%Usage:
%   mShed = BSLshedMask(Wsheds,Cut)
%       Wsheds = watershed masks for the region
%       Cut    = Specifies the number of sheds removed
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.%
%
%% BSL watershed mask calc

PT = Wsheds.PET;
Shed = Wsheds.Shed;
voxVol = Wsheds.voxVol;
[pX pY pZ] = size(PT);

% build shed masks
ShedTmp = zeros(size(PT));
indNZ = find(Shed);
Shed(isnan(Shed)) = 0;
Shed(isinf(Shed)) = 0;
ShedTmp(indNZ) = Shed(indNZ);

% removes very low background (i.e. bkg outside of patient)
SUVcut = 0.2;
ShedTmp(ShedTmp < SUVcut) = 0;
PT(ShedTmp <= 0) = 0;

% Shed trimming sucessively removing the lowest remaining shed
for k = 0:Cut
    minShed = min(nonzeros(ShedTmp(:)));
    indS = find(ShedTmp == minShed);
    ShedTmp(indS) = 0;
    PT(indS) = 0;
    if (isempty(indS) == 1)
        k = k - 1;
        continue
    end
end

if (sum(ShedTmp(:)) <=0)
    mPT = -1;
    return
end

mShed = ShedTmp;
mShed(mShed > 0) = 1;
voxVec = nonzeros(mShed.*PT);

% preliminary histogram estimate given the shed cuts
nVox = numel(voxVec);
SUViqr = iqr(voxVec);
IVHBinWidth = 2 * SUViqr * nVox^(-1/3);

[binHistV, volHistV]  = doseHist(voxVec, voxVol*ones(size(voxVec)), IVHBinWidth);
nBins = numel(binHistV);

% preliminary BKG estimate using histogram peak
[Hmax, iBkg] = max(volHistV);
Hbkg = binHistV(iBkg);
if (iBkg > 1 && iBkg < numel(binHistV))
    Hbkg = ( volHistV(iBkg-1)*binHistV(iBkg-1) ...
        + volHistV(iBkg)*binHistV(iBkg) ...
        + volHistV(iBkg+1)*binHistV(iBkg+1) ) ...
        / ( volHistV(iBkg-1) + volHistV(iBkg) + volHistV(iBkg+1) );
else
    if (iBkg < numel(binHistV))
        Hbkg = ( volHistV(iBkg)*binHistV(iBkg) ...
            + volHistV(iBkg+1)*binHistV(iBkg+1) ) ...
            / ( volHistV(iBkg) + volHistV(iBkg+1) );
    else
        Hbkg = ( volHistV(iBkg-1)*binHistV(iBkg-1) ...
            + volHistV(iBkg)*binHistV(iBkg) ) ...
            / ( volHistV(iBkg-1) + volHistV(iBkg) );
    end
end
% removes some activity from the left side of the histogram to remove
% fitting bias from activity below the background
SUVcut = Hbkg/3;
mShed(PT < SUVcut) = 0;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function printModelStats(model, printModelIssues, printDetails)
% printModelStats
%   prints some statistics about a model to the screen
%
%   model               a model structure
%   printModelIssues    true if information about unconnected
%                       reactions/metabolites and elemental balancing
%                       should be printed (opt, default false)
%   printDetails        true if detailed information should be printed
%                       about model issues. Only used if printModelIssues
%                       is true (opt, default true)
%
%   Usage: printModelStats(model,printModelIssues, printDetails)

if nargin<2
    printModelIssues=false;
end
if nargin<3
    printDetails=true;
end

fprintf(['Network statistics for ' model.id ': ' model.name '\n']);

%Get which reactions are present in each compartment
rxnComps=sparse(numel(model.rxns),numel(model.comps));

%For each compartment, find the metabolites that are present in that
%compartment and then the reactions they are involved in
for i=1:numel(model.comps)
    [~, I]=find(model.S(model.metComps==i,:));
    rxnComps(I,i)=1;
end

if isfield(model,'eccodes')
    fprintf(['EC-numbers\t\t\t' num2str(numel(unique(model.eccodes))) '\n']);
end

%Print information about genes
if isfield(model,'genes')
    fprintf(['Genes*\t\t\t\t' num2str(numel(model.genes)) '\n']);
    %Find the genes in each compartment
    for i=1:numel(model.comps)
        [~, I]=find(model.rxnGeneMat(rxnComps(:,i)==1,:));
        fprintf(['\t' model.compNames{i} '\t' num2str(numel(unique(I))) '\n']);
    end
end

%Print information about reactions
fprintf(['\nReactions*\t\t\t' num2str(numel(model.rxns)) '\n']);
for i=1:numel(model.comps)
    fprintf(['\t' model.compNames{i} '\t' num2str(sum(rxnComps(:,i))) '\n']);
end

%Removes the effect of compartments and removes duplicate reactions
temp=model;
temp.comps(:)={'s'}; %Set all compartments to be the same
equ=constructEquations(sortModel(temp,true,true),temp.rxns,false);

fprintf(['Unique reactions**\t' num2str(numel(unique(equ))) '\n']);

%Print information about metabolites
fprintf(['\nMetabolites\t\t\t' num2str(numel(model.mets)) '\n']);
for i=1:numel(model.comps)
    fprintf(['\t' model.compNames{i} '\t' num2str(sum(model.metComps==i)) '\n']);
end

fprintf(['Unique metabolites\t' num2str(numel(unique(model.metNames))) '\n']);

fprintf('\n* Genes and reactions are counted for each compartment if any of the corresponding metabolites are in that compartment. The sum may therefore not add up to the total number.\n');
fprintf('** Unique reactions are defined as being biochemically unique (no compartmentalization)\n');

%Also print some potential problems if there are any
if printModelIssues==true
    fprintf(['\nShort model quality summary for ' model.id ': ' model.name '\n']);

    %Check that all the metabolites are being used
    involvedMat=model.S;
    involvedMat(involvedMat~=0)=1;
    usedMets=sum(involvedMat,2);
    notPresent=find(usedMets==0);
    if ~isempty(notPresent)
        errorText=['Non-used metabolites\t' num2str(numel(notPresent)) '\n'];
        if printDetails==true
            for i=1:numel(notPresent)
                errorText=[errorText '\t(' model.mets{notPresent(i)} ') ' model.metNames{notPresent(i)} '\n'];
            end
            errorText=[errorText '\n'];
        end
        fprintf(errorText);
    end

    %Check if there are empty reactions
    usedRxns=sum(involvedMat,1);
    notUsed=find(usedRxns==0);
    if ~isempty(notUsed)
        errorText=['Empty reactions\t' num2str(numel(notUsed)) '\n'];
        if printDetails==true
            for i=1:numel(notUsed)
                errorText=[errorText '\t' model.rxns{notUsed(i)} '\n'];
            end
            errorText=[errorText '\n'];
        end
        fprintf(errorText);
    end

    %Check if there are dead-end reactions/metabolites
    [~, deletedReactions, deletedMetabolites]=simplifyModel(model,true,false,false,true);

    if ~isempty(deletedReactions)
        errorText=['Dead-end reactions\t' num2str(numel(deletedReactions)) '\n'];
        if printDetails==true
            for i=1:numel(deletedReactions)
                errorText=[errorText '\t' deletedReactions{i} '\n'];
            end
            errorText=[errorText '\n'];
        end
        fprintf(errorText);
    end

    %Ignore non-used metabolites
    deletedMetabolites=setdiff(deletedMetabolites,model.mets(notPresent));
    %Must map to indexes in order to print names
    deletedMetabolites=find(ismember(model.mets,deletedMetabolites));
    if ~isempty(deletedMetabolites)
        errorText=['Dead-end metabolites\t' num2str(numel(deletedMetabolites)) '\n'];
        if printDetails==true
            for i=1:numel(deletedMetabolites)
                errorText=[errorText '\t(' model.mets{deletedMetabolites(i)} ') ' model.metNames{deletedMetabolites(i)} '\n'];
            end
            errorText=[errorText '\n'];
        end
        fprintf(errorText);
    end

    balanceStructure=getElementalBalance(model);

    notParsed=find(balanceStructure.balanceStatus<0);
    notBalanced=find(balanceStructure.balanceStatus==0);

    if ~isempty(notParsed)
        errorText=['Reactions which could not be elementally balanced\t' num2str(numel(notParsed)) '\n'];
        if printDetails==true
            for i=1:numel(notParsed)
                errorText=[errorText '\t' model.rxns{notParsed(i)} '\n'];
            end
            errorText=[errorText '\n'];
        end
        fprintf(errorText);
    end
    if ~isempty(notBalanced)
        errorText=['Reactions which are elementally unbalanced\t' num2str(numel(notBalanced)) '\n'];
        if printDetails==true
            names=strcat(balanceStructure.elements.names,{', '});
            for i=1:numel(notBalanced)
                badOnes=sprintf('%s', names{abs(balanceStructure.leftComp(notBalanced(i),:)-balanceStructure.rightComp(notBalanced(i),:))>10^-7});
                errorText=[errorText '\t' model.rxns{notBalanced(i)} '\t' badOnes(1:end-2) '\n'];
            end
            errorText=[errorText '\n'];
        end
        fprintf(errorText);
    end
end
end

author: SysBioChalmers

repo: RAVEN


text:
function grf_io_test02 ( )

%*****************************************************************************80
%
%% GRF_IO_TEST02 tests GRF_HEADER_READ and GRF_DATA_READ.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 January 2009
%
%  Author:
%
%    John Burkardt
%
  input_filename = 'coxeter.grf';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'GRF_IO_TEST02\n' );
  fprintf ( 1, '  GRF_HEADER_READ reads the header of a GRF file.\n' );
  fprintf ( 1, '  GRF_DATA_READ reads the data of a GRF file.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Reading the GRF file "%s".\n', input_filename );

  [ node_num, edge_num ] = grf_header_read ( input_filename );

  grf_header_print ( node_num, edge_num );

  [ edge_pointer, edge_data, xy ] = grf_data_read ( input_filename, ...
    node_num, edge_num );

  grf_data_print ( node_num, edge_num, edge_pointer, edge_data, xy );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [irf_estimates,D,gamma,D_estimates,gamma_estimates,strshocks_estimates,favar]...
    =olsirft5(betahat,IRFperiods,Y,X,n,m,p,k,endo,pref,IRFband,names,enddate,startdate,T,data_endo,data_exo,const,strctident,IRFt,IRF,favar)

% computes IRF values (point estimates and confidence bands) for the OLS VAR model
% inputs:  - vector 'betahat': OLS VAR coefficients in vectorised form (defined in 1.1.15)
%          - matrix 'sigmahat': OLS VAR variance-covariance matrix of residuals (defined in 1.1.10)
%          - integer 'IRFperiods': number of periods for IRFs
%          - integer 'IRFt': determines which type of structural decomposition to apply (none, Choleski, triangular factorization)
%          - matrix 'X': matrix of regressors for the VAR model (defined in 1.1.8)
%          - integer 'n': number of endogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'm': number of exogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'p': number of lags included in the model (defined p 7 of technical guide)
%          - integer 'k': number of coefficients to estimate for each equation in the BVAR model (defined p 7 of technical guide)
%          - integer 'q': total number of coefficients to estimate for the BVAR model (defined p 7 of technical guide)
%          - scalar 'IRFband': confidence level for IRFs
% outputs: - cell 'irf_estimates': lower bound, point estimates, and upper bound for the IRFs
%          - matrix 'D': structural matrix for the OLS model (defined in 2.3.3)
%          - matrix 'gamma': structural disturbance variance-covariance matrix (defined p 48 of technical guide)
%          - vector 'D_estimates': point estimate (median) of the structural matrix D, in vectorised form
%          - vector 'gamma_estimates': point estimate (median) of the structural disturbance variance-covariance matrix gamma, in vectorised form

It=1001;
Bu=0;
Acc=It-Bu;

%% IV routine
[~,~,~,~,~,irf_storage,storage2]=...
    bear.irfIVols(betahat,n,m,p,k,T,names,startdate,enddate,X,Y,endo,data_endo,data_exo,const,pref,strctident,IRFt,IRFperiods,It,Bu);

%% reorganise storage
% loop over iterations
for ii=1:Acc
    % loop over IRF periods
    for jj=1:IRFperiods
        % loop over variables
        for kk=1:n
            % loop over shocks
            for ll=1:n
                irf_record{kk,ll}(ii,jj)=irf_storage{ii,1}(kk,ll,jj);
            end
        end
    end
    D_record(:,ii)=storage2{ii,1}(:);
end

for ii=1:n^2
    D_estimates(ii,1)=quantile(D_record(ii,:),0.5);
end
D=reshape(D_estimates,n,n);
D_estimates=D(:);
gamma=eye(n); %equivalent to the gamma_record output of irfres
gamma_estimates=gamma(:);

% create then the cell storing the point estimates and confidence bands
irf_estimates=cell(n,n);

if IRF==1 | favar.IRFplot==1
    % rearrange
    for ii=1:n
        % loop over variables
        for jj=1:n
            % loop over IRF periods
            for kk=1:IRFperiods
                % median
                irf_estimates{jj,ii}(2,kk)=quantile(irf_record{jj,ii}(:,kk),0.5);
                % lower bound
                irf_estimates{jj,ii}(1,kk)=quantile(irf_record{jj,ii}(:,kk),(1-IRFband)/2);
                % upper bound
                irf_estimates{jj,ii}(3,kk)=quantile(irf_record{jj,ii}(:,kk),IRFband+(1-IRFband)/2);
            end
        end
    end
end

% check if the variables have been transformed
if favar.FAVAR==1
    if favar.transformation==1 || favar.plot_transform==1
        % re-transform irf_estimates
        favar.IRF.irf_estimates_nottransformed=irf_estimates; % before, save untransformed IRFs
        % re-transform
        [irf_estimates]=bear.favar_retransX_irf_estimates(irf_estimates,favar.transformationindex_endo,favar.levels);
    end
end

% finally, estimate structural shocks
% Compute first the model residuals
Bhat=reshape(betahat,k,n);
EPS=Y-X*Bhat;
% Then use (XXX) to recover the structural shocks
ETA=D\EPS';
% output
strshocks_estimates=ETA;

author: european-central-bank

repo: BEAR-toolbox


text:
function [coef, bias, variance] = regressMajorityPercentage(smaps, pv, ph, imsegs)
% estimate the percentage of a region occupied by the majority label, given
% some estimate of the superpixel label likelihoods

ndata = 0;
for f = 1:numel(smaps)
    for m = 1:size(smaps{f}, 2)
        ndata = ndata + max(smaps{f}(:, m));
    end
    plab{f} = [pv{f}(:, 1) repmat(pv{f}(:, 2), [1 5]).*ph{f} pv{f}(:, 3)];
end

disp(num2str(ndata))

y = zeros(ndata, 1); % label
x = ones(ndata, 3); % statistics for regression
c = 0;
if 1
for f = 1:numel(smaps)
    npix = imsegs(f).npixels(:);
    labels = imsegs(f).labels(:);

    for m = 1:size(smaps{f}, 2)
        for s = 1:max(smaps{f}(:, m))
            ind = find(smaps{f}(:, m)==s);

            mclab = 0;
            mcprc = 0;
            for k = 1:7
                prc = sum((labels(ind)==k).*npix(ind))/(sum(npix(ind).*(labels(ind)>0))+1E-10);
                if prc > mcprc
                    mcprc = prc;
                    mclab = k;
                end
            end
            if mclab > 0
                y(c+s) = mcprc;
                x(c+s, 1) = sum(plab{f}(ind, mclab).*npix(ind))/sum(npix(ind));
                [tmp, maxlab] = max(plab{f}(ind, :), [], 2);
                x(c+s, 2) = sum((maxlab==mclab).*npix(ind)) / sum(npix(ind));
                x(c+s, 3) = numel(ind);
            end
        end
        c = c + s;
    end
end
ind = find(y==0);
y(ind) = [];
x(ind, :) = [];
save '../data/tmp.mat' x y
else
    load '../data/tmp.mat'
end

thresh1 = [0.05:0.1:0.95];
tx = sum(repmat(x(:, 1), [1 numel(thresh1)])  >= repmat(thresh1, [size(x,1) 1]), 2);
hist(tx)
for k =1:numel(thresh1)
    ind = find(tx==k);
    px(k) = mean(y(ind));
    count(k) = numel(ind);
    x(ind, 1) = px(k);
end
disp(num2str(px))
disp(num2str(count))

thresh2 = [0.05:0.1:0.95];
tx = sum(repmat(x(:, 2), [1 numel(thresh2)])  >= repmat(thresh2, [size(x,1) 1]), 2);
hist(tx)
for k =1:numel(thresh2)
    ind = find(tx==k);
    px(k) = mean(y(ind));
    count(k) = numel(ind);
    x(ind, 2) = px(k);
end
disp(num2str(px))
disp(num2str(count))

thresh3 = [1 2 3 4 5 7 10 15 20 25 35 50 100 150];
tx = sum(repmat(x(:, 3), [1 numel(thresh3)])  >= repmat(thresh3, [size(x,1) 1]), 2);
hist(tx)
for k =1:numel(thresh3)
    ind = find(tx==k);
    px(k) = mean(y(ind));
    count(k) = numel(ind);
    x(ind, 3) = px(k);
end
disp(num2str(px))
disp(num2str(count))

ylog = log((y+1E-5)./(1-y+1E-5));
xlog = log((x+1E-5)./(1-x+1E-5));

coef = fminunc(@(a) objective(a, xlog, y), [1 1 1 0]');

%coef = robustfit(xlog, ylog, 'logistic');
%coef = regress(ylog, [ones([size(xlog, 1) 1])  xlog]);
%coef = regress(ylog, xlog);

%ylogEst = coef(1) + xlog*coef(2:end);
ylogEst = coef(end) + xlog*coef(1:end-1);
size(ylogEst)
yEst = 1./(1+exp(-ylogEst));
disp(num2str([mean(yEst) mean(y)]))
bias = mean(yEst-y);
variance = var(yEst-y);
disp(mean(abs(yEst-y)))

function val = objective(coef, xlog, y)
ylogEst = coef(end) + xlog*coef(1:end-1);
yEst = 1./(1+exp(-ylogEst));
val = sum((yEst-y).^2);
%disp([sqrt(val/numel(y)) coef'])

author: zouchuhang

repo: LayoutNet


text:
function [dfdz,dfdp,sym_dfdz,sym_dfdp] = symDynJac(ode,nstates,nparam)
%SYMJAC  Returns Symbolically Differentiated Partial Derivatives of a Dynamic System
%
%   [dfdz,dfdp] = symDynJac(ode) uses the symbolic toolbox to automatically
%   generate the sensitivity partial derivatives of the function handle ode.
%   You should supply the function handle in the form @(t,z,p) z(1)^2 + p(1)
%   noting t, z and p must be the only variables used, and are indexed in the
%   equation (no vector operations).
%
%   [dfdz,dfdp] = symDynJac(ode,nstates) specifies the number of states in
%   the equation, assuming consecutive ordering. Useful if a state is not
%   specified in the original equation to pad DFDZ with zeros.
%
%   [dfdz,dfdp] = symDynJac(ode,nstates,nparam) specifies the number of
%   parameters in the equation, assuming consecutive ordering. Useful if a
%   parameter is not specified in the original equation to pad DFDP with
%   zeros.

%   Copyright (C) 2013 Jonathan Currie (IPL)

if(nargin < 3), nparam = 0; end
if(nargin < 2), nstates = 0; end

if(~optiCheckSymTBX())
    dfdz = []; sym_dfdz = [];
    dfdp = []; sym_dfdp = [];
    return
end

if(~isa(ode,'function_handle') && ~isa(ode,'barvec'))
    error('Fun should be a function handle!');
end
if(isa(ode,'function_handle') && nargin(ode) ~= 3)
    error('ODE should only have three input arguments (t,z,p)');
end

%Convert ODE to a symbolic expression
if(isa(ode,'function_handle'))
    [symode,ind] = func2sym(ode,{'z','p','t'});
    indz = ind{1}; indp = ind{2};
else
    symode = sym(getEq(ode)); %convert from barvec to symbolic expression
    v = char(symvar(symode));
    indz = 1:length(strfind(v,'z'));
    indp = 1:length(strfind(v,'p'));
end

%Create each symbolic partial derivative
sym_dfdz = symPartialDer(symode,'z',nstates,indz);
sym_dfdp = symPartialDer(symode,'p',nparam,indp);

%Return to function handles (note order of args important and changed)
dfdz = sym2func(sym_dfdz,{'t','z','p'});
dfdp = sym2func(sym_dfdp,{'t','z','p'});

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% Test file for bndfun/restrict.m

function pass = test_restrict(pref)

% Get preferences.
if (nargin < 1)
    pref = chebfunpref();
end

% Set the domain
dom = [-2 7];

%%
% Check behavior for empty inputs.
f = bndfun();
f = restrict(f, [-0.5 0.5]);
pass(1) = isempty(f);

%%
% Check behaviour for non-subinterval inputs.
f = bndfun(@(x) sin(x), struct('domain', dom), pref);
g = restrict(f, dom);
pass(2) = isequal(f, g);

try
    g = restrict(f, [dom(1) - 1, dom(2) + 1]); %#ok<NASGU>
    pass(3) = 0;
catch ME
    pass(3) = strcmp(ME.identifier, 'CHEBFUN:BNDFUN:restrict:badInterval');
end

try
    g = restrict(f, [-Inf, 1]); %#ok<NASGU>
    pass(4) = 0;
catch ME
    pass(4) = strcmp(ME.identifier, 'CHEBFUN:BNDFUN:restrict:badInterval');
end

try
    g = restrict(f, [-1 -0.25 0.3 0.1 1]); %#ok<NASGU>
    pass(5) = 0;
catch ME
    pass(5) = strcmp(ME.identifier, 'CHEBFUN:BNDFUN:restrict:badInterval');
end

%%
% Check whether restriction actually results in a BNDFUN on the correct domain.
g = restrict(f,[2, 3]);
pass(6) = all(g.domain == [2, 3]);

%%
% Spot-check a few functions
pass(7) = test_spotcheck_restrict(@(x) exp(x) - 1, dom, [-2 4], [], pref);
pass(8) = test_spotcheck_restrict(@(x) 1./(1 + x.^2), dom, [-0.7 0.9], [], ...
    pref);
pass(9) = test_spotcheck_restrict(@(x) cos(1e3*x), dom, [0.1 0.5], [], pref);
pass(10) = test_spotcheck_restrict(@(t) sinh(t*exp(2*pi*1i/6)), dom, ...
    [-0.4 1], [], pref);

%%
% Check multiple subinterval restriction.
f = bndfun(@(x) sin(x) + sin(x.^2), struct('domain', dom), pref);
g = restrict(f, [-1.7 2.3 6.8]);
h1 = restrict(f, [-1.7 2.3]);
h2 = restrict(f, [2.3 6.8]);
x = linspace(-1, 1, 100).';
err1 = norm(feval(g{1} - h1, x), inf);
err2 = norm(feval(g{2} - h2, x+4), inf);
tol = 1e3*eps;

pass(11) = err1 < tol && err2 < tol;

%%
% Check whether restriction actually results in a BNDFUN on the correct domain.
g = restrict(f,[2, 3, 5]);
pass(12) = all(g{1}.domain == [2, 3] & g{2}.domain == [3,5]);

%%
% Check operation for array-valued functions.
pass(13) = test_spotcheck_restrict(@(x) [sin(x) cos(x) exp(x)], dom, ...
    [-1 -0.7], [], pref);

f = bndfun(@(x) [sin(x) cos(x)], struct('domain', dom), pref);
g = restrict(f, [-0.6 0.1 1]);
h1 = restrict(f, [-0.6 0.1]);
h2 = restrict(f, [0.1 1]);
x = linspace(-.5, 0, 100).';
err1 = norm(feval(g{1} - h1, x), inf);
err2 = norm(feval(g{2} - h2, x+.7), inf);
tol = 10*eps;
pass(14) = err1 < tol && err2 < tol;

%% Test on singular function:

pow = -0.5;
op = @(x) (x - dom(1)).^pow.*sin(x);
pref.blowup = true;
data.exponents = [pow 0];
pass(15) = test_spotcheck_restrict(op, dom, [-1 -0.7], data, pref);

end

% Spot-check restriction of a given function to a given subinterval.
function result = test_spotcheck_restrict(fun_op, dom, subint, data, pref)
% Perform restriction.
data.domain = dom;
f = bndfun(fun_op, data, pref);
g = restrict(f, subint);

a = subint(1);
b = subint(2);

% Sample on a grid of 100 points and check for accuracy.
x = linspace(a, b, 100).';
y_exact = fun_op(x);
y_approx = feval(g, x);

result = norm(y_exact - y_approx, Inf) <...
         8e4*max(get(f, 'vscale')*eps);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function m = size(t,idx)
%SIZE Sparse tensor dimensions.
%
%   D = SIZE(T) returns the size of the tensor.
%
%   I = size(T,DIM) returns the sizes of the dimensions specified by DIM.
%
%   See also SPTENSOR, SPTENSOR/NDIMS.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

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

if exist('idx','var')
    m = t.size(idx);
else
    m = t.size;
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
%% Main function to generate tests
function tests = testSuite

tests = functiontests(localfunctions);

end

%% Test Functions
function wgs2gcjTest(testCase)
% test wgs to gcj

TESTS =  [
    31.1774276, 121.5272106, 31.17530398364597, 121.531541859215; %shanghai
    22.543847, 113.912316, 22.540796131694766, 113.9171764808363; %shenzhen
    39.911954, 116.377817, 39.91334545536069, 116.38404722455657; %beijing
    39.739200,-104.990300,39.739200,-104.990300; %Denver edge case
    -27.469800,153.025100, -27.469800,153.025100; % Brisbane edge case
    -22.906800,-43.172900, -22.906800,-43.172900; % Rio de Janerio edge case
    ];

[actLat actLong] = wgs2gcj(TESTS(:,1),TESTS(:,2));
expLat = TESTS(:,3);
expLong = TESTS(:,4);

verifyEqual(testCase,actLat,expLat, 'AbsTol' , .000001);
verifyEqual(testCase,actLong,expLong, 'AbsTol' , .000001);

end

function gcj2wgsTest(testCase)
% test wgs to gcj

TESTS =  [
    31.1774276, 121.5272106, 31.17530398364597, 121.531541859215; %shanghai
    22.543847, 113.912316, 22.540796131694766, 113.9171764808363; %shenzhen
    39.911954, 116.377817, 39.91334545536069, 116.38404722455657; %beijing
    39.739200,-104.990300,39.739200,-104.990300; %Denver edge case
    -27.469800,153.025100, -27.469800,153.025100; % Brisbane edge case
    -22.906800,-43.172900, -22.906800,-43.172900; % Rio de Janerio edge case
    ];

[actLat actLong] = gcj2wgs(TESTS(:,3),TESTS(:,4));
expLat = TESTS(:,1);
expLong = TESTS(:,2);

measurementDistance = distance(expLat, expLong, actLat,actLong);
actualDistance = zeros(length(measurementDistance),1);
verifyEqual(testCase,actualDistance,measurementDistance, 'AbsTol' , 5);

end

function gcj2wgsExactTest(testCase)
% test wgs to gcj

TESTS =  [
    31.1774276, 121.5272106, 31.17530398364597, 121.531541859215; %shanghai
    22.543847, 113.912316, 22.540796131694766, 113.9171764808363; %shenzhen
    39.911954, 116.377817, 39.91334545536069, 116.38404722455657; %beijing
    39.739200,-104.990300,39.739200,-104.990300; %Denver edge case
    -27.469800,153.025100, -27.469800,153.025100; % Brisbane edge case
    -22.906800,-43.172900, -22.906800,-43.172900; % Rio de Janerio edge case
    ];

[actLat actLong] = gcj2wgs_exact(TESTS(:,3),TESTS(:,4));
expLat = TESTS(:,1);
expLong = TESTS(:,2);

measurementDistance = distance(expLat, expLong, actLat,actLong);
actualDistance = zeros(length(measurementDistance),1);
verifyEqual(testCase,actualDistance,measurementDistance, 'AbsTol' , .5);

end

function mars2wgsTest(testCase)
% test wgs to gcj should be very exact and very fast....

TESTS =  [
    31.1774276, 121.5272106, 31.17530398364597, 121.531541859215; %shanghai
    22.543847, 113.912316, 22.540796131694766, 113.9171764808363; %shenzhen
    39.911954, 116.377817, 39.91334545536069, 116.38404722455657; %beijing
    39.739200,-104.990300,39.739200,-104.990300; %Denver edge case
    -27.469800,153.025100, -27.469800,153.025100; % Brisbane edge case
    -22.906800,-43.172900, -22.906800,-43.172900; % Rio de Janerio edge case
    ];

[actLat actLong] = mars2wgs(TESTS(:,3),TESTS(:,4));
expLat = TESTS(:,1);
expLong = TESTS(:,2);

measurementDistance = distance(expLat, expLong, actLat,actLong);
actualDistance = zeros(length(measurementDistance),1);
verifyEqual(testCase,actualDistance,measurementDistance, 'AbsTol' , .01);

end

function wgs2bdTest(testCase)
% test wgs to gcj should be very exact and very fast....

 TESTS_BD = [
    29.199786, 120.019809, 29.196131605295484, 120.00877901149691;
    29.210504, 120.036455, 29.206795749156136, 120.0253853970846
    ];

[actLat actLong] = wgs2bd(TESTS_BD(:,3),TESTS_BD(:,4));
expLat = TESTS_BD(:,1);
expLong = TESTS_BD(:,2);

verifyEqual(testCase,actLat,expLat, 'AbsTol' , .00005);
verifyEqual(testCase,actLong,expLong, 'AbsTol' , .00005);

end

function bd2wgsTest(testCase)
% test wgs to gcj should be very exact and very fast....

 TESTS_BD = [
    29.199786, 120.019809, 29.196131605295484, 120.00877901149691;
    29.210504, 120.036455, 29.206795749156136, 120.0253853970846
    ];

[actLat actLong] = bd2wgs(TESTS_BD(:,1),TESTS_BD(:,2));
expLat = TESTS_BD(:,3);
expLong = TESTS_BD(:,4);

verifyEqual(testCase,actLat,expLat, 'AbsTol' , .00005);
verifyEqual(testCase,actLong,expLong, 'AbsTol' , .00005);

end

function gcj2bdTest(testCase)
% test wgs to gcj should be very exact and very fast....
% bd lat, bd long, gcj lat, gcj long
 TESTS_BD = [39.90851,116.43351,39.90245,116.42703];

[actLat actLong] = gcj2bd(TESTS_BD(:,3),TESTS_BD(:,4));
expLat = TESTS_BD(:,1);
expLong = TESTS_BD(:,2);

verifyEqual(testCase,actLat,expLat, 'AbsTol' , .00005);
verifyEqual(testCase,actLong,expLong, 'AbsTol' , .00005);

end

function bd2gcjTest(testCase)
% test wgs to gcj should be very exact and very fast....
% bd lat, bd long, gcj lat, gcj long
 TESTS_BD = [39.90851,116.43351,39.90245,116.42703];

[actLat actLong] = bd2gcj(TESTS_BD(:,1),TESTS_BD(:,2));
expLat = TESTS_BD(:,3);
expLong = TESTS_BD(:,4);

verifyEqual(testCase,actLat,expLat, 'AbsTol' , .00005);
verifyEqual(testCase,actLong,expLong, 'AbsTol' , .00005);

end

function deg2utmTest(testCase)
% test wgs to utm should be very exact and very fast....
% LAT, LONG, X, Y

TestUTMDEG = [36.903784,-104.545898,540455.18,4084295.17];

[actUTMx, actUTMy] = deg2utm(TestUTMDEG(:,1),TestUTMDEG(:,2));
expUTMx = TestUTMDEG(:,3);
expUTMy = TestUTMDEG(:,4);

verifyEqual(testCase,actUTMx,expUTMx, 'AbsTol' , .1);
verifyEqual(testCase,actUTMy,expUTMy, 'AbsTol' , .1);

end

function utm2degTest(testCase)
% test wgs to gcj should be very exact and very fast....
% LAT, LONG, X, Y

TestUTMDEG = [36.903784,-104.545898,540455.18,4084295.17];

[actLat, actLong] = utm2deg(TestUTMDEG(:,3),TestUTMDEG(:,4),'13N');
expLat = TestUTMDEG(:,1);
expLong = TestUTMDEG(:,2);

verifyEqual(testCase,actLat,expLat, 'AbsTol' , .00001);
verifyEqual(testCase,actLong,expLong, 'AbsTol' , .00001);

end

function speedTest(testCase)
% test wgs to gcj
n = 50000;
TESTS =  [
    31.1774276, 121.5272106, 31.17530398364597, 121.531541859215; %shanghai
    22.543847, 113.912316, 22.540796131694766, 113.9171764808363; %shenzhen
    39.911954, 116.377817, 39.91334545536069, 116.38404722455657; %beijing
    39.739200,-104.990300,39.739200,-104.990300; %Denver edge case
    -27.469800,153.025100, -27.469800,153.025100; % Brisbane edge case
    -22.906800,-43.172900, -22.906800,-43.172900; % Rio de Janerio edge case
    ];

tic;
for i=1:n
    [actLat actLong] = wgs2gcj(TESTS(:,1),TESTS(:,2));
end
times.wgs2gcj = toc;

tic;
for i=1:n
    [actLat actLong] = mars2wgs(TESTS(:,1),TESTS(:,2));
end
times.mars2gcj = toc;

tic;
for i=1:n
    [actLat actLong] = gcj2wgs(TESTS(:,1),TESTS(:,2));
end
times.gcj2wgs = toc;

tic;
for i=1:n
    [actLat actLong] = gcj2wgs_exact(TESTS(:,1),TESTS(:,2));
end
times.gcj2wgs_exact = toc;

times

end

 function setupOnce(testCase)  % do not change function name

    testsBD = [
    29.199786, 120.019809, 29.196131605295484, 120.00877901149691;
    29.210504, 120.036455, 29.206795749156136, 120.0253853970846
    ];

 end

 function teardownOnce(testCase)  % do not change function name
% change back to original path, for example
end

function setup(testCase)  % do not change function name
% open a figure, for example
end

function teardown(testCase)  % do not change function name
% close figure, for example
end

author: googollee

repo: eviltransform


text:
function chebyshev_polynomial_test09 ( )

%*****************************************************************************80
%
%% CHEBYSHEV_POLYNOMIAL_TEST09 compares a function and projection over [-1,+1].
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 April 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CHEBYSHEV_POLYNOMIAL_TEST09:\n' );
  fprintf ( 1, '  T_PROJECT_COEFFICIENTS computes the Chebyshev interpolant C(F)(N,X)\n' );
  fprintf ( 1, '  of a function F(X) defined over [-1,+1].\n' );
  fprintf ( 1, '  T_PROJECT_VALUE evaluates that projection.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Compute projections of order N to exp(x) over [-1,+1],\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   N   Max||F(X)-C(F)(N,X)||\n' );
  fprintf ( 1, '\n' );
  for n = 0 : 10
    c = t_project_coefficients ( n, @exp );
    m = 101;
    x = ( linspace ( -1.0, +1.0, m ) )';
    v = t_project_value ( m, n, x, c );
    r = v - exp ( x );
    fprintf ( 1, '  %2d  %12.4g\n', n, max ( abs ( r ) ) );
  end

  return
end


text:
function [mNewCatalog] = bc_SelectMechanism_Start(mCatalog)

% Default return value
mNewCatalog = nan;

% Invoke the user interface
hDialog = bc_SelectMechanism_Dialog;

% Analyze Output
if ~ishandle(hDialog)
  answer = 0;
else
  handles = guidata(hDialog);
  answer = handles.answer;
  % OK pressed
  if answer == 1
    % Get the user options
    nType = get(handles.cboMechanism, 'Value');
    fAngle = str2double(get(handles.txtAngle, 'String'));
    % Remove figure from memory
    delete(hDialog);
    % Create the new catalog
    [mNewCatalog] = bc_SelectMechanism(mCatalog, nType, fAngle);
  else
    % Remove figure from memory
    delete(hDialog);
  end
end


text:
function [ranset,seed]=randomset(tot,nel,seed)

%RANDOMSET draws randomly a subsample of nel cases out of tot.
%(It is called if not all (p+1)-subsets out of n will be considered.)
%
% Required input arguments:
%     tot : The total number of observations to consider
%     nel : The number of observations that the subsample must contain
% Optional input arguments:
%     seed : To define the state of the generator (default=0)
%            (0 sets the generator to its default initial state)
%
% Output arguments:
%   ranset : Random subset of nel cases out of tot.
%   seed   : The corresponding state.
%
%
% I/O:
%    [ranset,seed]=randomset(n,n/2,0)
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at:
%              http://wis.kuleuven.be/stat/robust.html

if nargin==2
    seed=0;
end

for j=1:nel
   [random,seed]=uniran(seed);
   num=floor(random*tot)+1;
   if j > 1
      while any(ranset==num)
         [random,seed]=uniran(seed);
         num=floor(random*tot)+1;
      end
   end
   ranset(j)=num;
end

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [fx,dF_dX,dF_dTheta] = f_doubleWell(Xt,Theta,ut,inF)
% damped double Well Lyapunov evolution function
% function [f,J] = f_doubleWell(t,x,theta)
%
% This function computes the evolution function that comes from a .
% IN:
%   - t: time index (not used here)
%   - x: the current state of the system
%   - theta: a 2x1 vector parameter (containing the position of the two
%   wells.
% OUT:
%   - f: the current value of the evolution function
%   - J: the jacobbian of the system

deltat = inF.deltat;

a   = Theta(1);
b   = Theta(2);
k   = Theta(3);
x1  = Xt(1);
x2  = Xt(2);

f   = [ x2  ;  -2*(x1-a)*(x1-b)^2-2*(x1-a)^2*(x1-b)-k.*x2 ];
J   = [ 0                                               1
        -2*(x1-b)^2-8*(x1-a)*(x1-b)-2*(x1-a)^2        -k ];

fx                  = deltat.*f + Xt;
dF_dX               = deltat.*J' + eye(2);
dF_dTheta           = deltat.*[ 0 2*(x1-b)^2+4*(x1-a)*(x1-b)
                                0 4*(x1-a)*(x1-b)+2*(x1-a)^2
                                0 -x2                        ];

author: MBB-team

repo: VBA-toolbox


text:
function gmt_scale(hAxes, sOutput, sPrefix)

% Get the figure-handle
hFigure = get(hAxes, 'Parent');

% Is there a colorbar?
vKids = get(hFigure, 'Children');
for nCnt = length(vKids):-1:1
  sType = get(vKids(nCnt), 'Type');
  sTag = get(vKids(nCnt), 'Tag');
  if (strcmp(sType, 'axes') & strcmp(sTag, 'Colorbar'))
    vUserData = get(vKids(nCnt), 'UserData');
    if vUserData.PlotHandle == hAxes

    % Extract the colormap
     vCLim = get(hAxes, 'CLim');
     mColormap = colormap(hAxes);
     [nRow, nCol] = size(mColormap);
     fRange = vCLim(2)-vCLim(1);
     fStep = fRange/(nRow-1);

     % Create the cpt-file and save it
     mMap = [];
     for nRowCnt = 1:nRow-1
       mMap = [mMap; (vCLim(1)+((nRowCnt-1)*fStep)) mColormap(nRowCnt,:).*255 (vCLim(1)+(nRowCnt*fStep)) mColormap(nRowCnt+1,:).*255];
     end
     save([sPrefix 'colormap.cpt'], 'mMap', '-ascii');

    % Get properties of colorbar
    vXTick = get(vKids(nCnt), 'XTick');
    vYTick = get(vKids(nCnt), 'YTick');
    if isempty(vXTick)
      % Colorbar is vertical
      sD = '';
  vYDiff = diff(vYTick);
 vSel = ~(vYDiff == 0);
 vYDiff = vYDiff(:,vSel);
 fYDiff = min(abs(vYDiff));
     sB = num2str(fYDiff);
    else
      % Colorbar is horizontal
      sD = 'h';
     vXDiff = diff(vXTick);
    vSel = ~(vXDiff == 0);
    vXDiff = vXDiff(:,vSel);
    fXDiff = min(abs(vXDiff));
     sB = num2str(fXDiff);
    end
%

    %get(vKids(nCnt))

    % Plot the colorbar
hFile = fopen(sOutput, 'a');
    fprintf(hFile, 'psscale -D6.3i/1.35i/3i/0.3i%s -C%s -B%s -O -K >> $output\n', sD, [sPrefix 'colormap.cpt'], sB);

fclose(hFile);

end
  end
end


text:
%rf_freqshift.m
%Jamie Near, McGill University 2014.
%
% USAGE:
% RF_shift=rf_freqshift(RF,Tp,F);
%
% DESCRIPTION:
% Apply a frequency shift to an RF pulse.
%
% INPUTS:
% RF         = RF pulse definition structure.
% Tp         = duration of the rf pulse in [ms].
% F          = amount that you would like to frequency shift the rf pulse in [Hz].
%
% OUTPUTS:
% RF_shift   = Output rf pulse following frequency shift.

function RF_shift=rf_freqshift(RF,Tp,F)

N=size(RF.waveform,1);
Tp=Tp/1000;
dt=Tp/N;
t=[0:dt:Tp-dt];

phaseRamp = t * F * 360;

RF_shift=RF;
RF_shift.waveform(:,1)=RF_shift.waveform(:,1)+phaseRamp';
RF_shift.f0=RF.f0+F;

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function [slug] = kg2slug(kg)
% Convert mass from kilograms to slugs.
% Chad Greene 2012
slug = kg/14.5939;

author: Sable

repo: mcbench-benchmarks


text:
function [NF,NL,NQc,NH,NPinf,NdF,NdQc,NdPinf,T] = ss_balance(F,L,Qc,H,Pinf,dF,dQc,dPinf)
% SS_BALANCE - Balance state space model for improved numerical stability
%
% Syntax:
%   [F,L,Qc,H,Pinf,dF,dQc,dPinf,T] = ss_balance(F,L,Qc,H,Pinf,dF,dQc,dPinf)
%
% In:
%   F           - Feedback matrix
%   L           - Noise effect matrix
%   Qc          - Spectral density of white noise process w(t)
%   H           - Observation model matrix
%   Pinf        - Covariance of the stationary process
%   dF          - Derivatives of F w.r.t. parameters
%   dQc         - Derivatives of Qc w.r.t. parameters
%   dPinf       - Derivatives of Pinf w.r.t. parameters
%
% Out:
%   (...)       - As above, but balanced
%   T           - T is the balancing matrix from 'balance'
%
% Description:
%   This function takes the SDE state space model matrices as inputs and
%   outputs the same matrices, but in a numerically more stable form.
%   This balancing is based on the Matlab 'balance' function (see the
%   reference).
%
%   The state space model is given as follows in terms of a stochastic
%   differential equation
%
%      df(t)/dt = F f(t) + L w(t),
%
%   where w(t) is a white noise process with spectral denisty Qc. The
%   observation model matrix is denoted by matrix H.
%
% References:
%   [1] Beresford N. Parlett and Christian Reinsch (1969). Balancing
%       a matrix for calculation of eigenvalues and eigenvectors.
%       Numerische Mathematik, 13(4): 293-304.
%
% See also:
%   BALANCE
%
% Copyright:
%   2014 Arno Solin and Simo Sarkka
%
%  This software is distributed under the GNU General Public
%  License (version 3 or later); please refer to the file
%  License.txt, included with the software, for details.

%% Balance the state space model

  % This is based on the following:
  %
  %  dx/dt = F x + L w
  %      y = H x
  %
  % Let T z = x, which gives
  %
  %  dz/dt = inv(T) F T z + inv(T) L w
  %      y = H T z
  %
  % See also: [A,B,C,E,S,P] = aebalance(A,B,C,E)

  % Balance the dynamic model matrix
  [T,NF] = balance(F);

  % Balance noise effect matrix
  NL = T\L;

  % Balance the measurement model
  NH = H*T;

  % Balance spectral density
  NQc = Qc;

  % Balance stationary state covariance matrix
  if nargin < 5
    NPinf = [];
  else
    L = chol(Pinf,'lower');
    LL = T\L;
    NPinf = LL*LL';
    %NPinf = T\Pinf/T;
  end

  % Balance partial derivatives in F
  if nargin < 6
    NdF = [];
  else
    NdF = dF;
    for j=1:size(dF,3)
      NdF(:,:,j) = T\dF(:,:,j)*T;
    end
  end

  % Balane partial derivatives in dQc;
  if nargin < 7
    NdQc = [];
  else
    NdQc = dQc;
  end

  % Balance partial derivatives in dF
  if nargin < 8
    NdPinf = [];
  else
    NdPinf = dPinf;
    for j=1:size(dF,3)
      NdPinf(:,:,j) = T\dPinf(:,:,j)/T;
    end
  end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function sVF = plus(sVF1, sVF2)
%
% Syntax
%   sVF = sVF1+sVF2
%   sVF = a+sVF1
%   sVF = sVF1+a
%

if isa(sVF1, 'vector3d')
  sVF = sVF2;
  sVF.x = sVF.x+sVF1.x;
  sVF.y = sVF.y+sVF1.y;
  sVF.z = sVF.z+sVF1.z;

elseif isa(sVF2, 'vector3d')
  sVF = sVF1;
  sVF.x = sVF.x+sVF2.x;
  sVF.y = sVF.y+sVF2.y;
  sVF.z = sVF.z+sVF2.z;

else
  sVF = sVF1;
  sVF.sF = sVF1.sF+sVF2.sF;

end

end

author: mtex-toolbox

repo: mtex

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


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

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

% Last Modified by GUIDE v2.5 08-Jun-2009 20:35:41

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

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

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

% Choose default command line output for About
handles.output = hObject;
% Update handles structure
movegui(hObject,'onscreen')% To display application onscreen
movegui(hObject,'center')  % To display application in the center of screen
% Update handles structure
guidata(hObject, handles);

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

function varargout = About_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

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

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

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

% Hint: delete(hObject) closes the figure
delete(hObject);

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

% --- Executes during object creation, after setting all properties.
function logo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to logo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called
a=imread('im3.png');
imshow(a);
% Hint: place code in OpeningFcn to populate logo


text:
function linpack_d_test02 ( )

%*****************************************************************************80
%
%% TEST02 tests DCHEX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2009
%
%  Author:
%
%    John Burkardt
%
  n = 5;
  lda = n;
  nz = 1;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST02\n' );
  fprintf ( 1, '  For double precision, general storage,\n' );
  fprintf ( 1, '  DCHEX can shift columns in a Cholesky factorization.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The number of equations is N = %d\n', n );
%
%  Set the values of the matrix A.
%
  a(1:n,1:n) = 0.0;

  for i = 1 : n
    a(i,i) = 2.0;
  end

  for i = 1 : n-1
    a(i,i+1) = -1.0;
  end

  for i = 2 : n
    a(i-1,i) = -1.0;
  end

  for i = 1 : n
    z(i,1) = i;
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The matrix A:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    for j = 1 : n
      fprintf ( 1, '  %14f', a(i,j) );
    end
    fprintf ( 1, '\n' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The vector Z:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    fprintf ( 1, '  %14f\n', z(i) );
  end
%
%  Decompose the matrix.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Decompose the matrix.\n' );

  job = 0;
  ipvt(1:n) = 0;

  [ a, ipvt, info ] = dchdc ( a, lda, n, ipvt, job );

  if ( info ~= n )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'DCHDC returned INFO = %d\n', info );
    fprintf ( 1, '  This means the matrix is not positive definite.\n' );
  end
%
%  Zero out the lower diagonal.
%
  for i = 2 : n
    for j = 1 : i-1
      a(i,j) = 0.0;
    end
  end
%
%  Print the factorization.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The Cholesky factor U:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    for j = 1 : n
      fprintf ( 1, '  %14f', a(i,j) );
    end
    fprintf ( 1, '\n' );
  end
%
%  Right circular shift columns L through K.
%
  k = 1;
  l = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, ...
    '  Right circular shift columns K  = %d through L = %d\n', k, l );

  job = 1;
  [ a, z, c, s ] = dchex ( a, lda, n, k, l, z, n, nz, job );
%
%  Left circular shift columns K+1 through L.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, ...
    '  Left circular shift columns K+1 = %d through L = %d\n', k+1, l );

  job = 2;
  [ a, z, c, s ] = dchex ( a, lda, n, k+1, l, z, n, nz, job );
%
%  Print the factorization.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The shifted Cholesky factor U:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    for j = 1 : n
      fprintf ( 1, '  %14f', a(i,j) );
    end
    fprintf ( 1, '\n' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The shifted vector Z:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    fprintf ( 1, '  %14f\n', z(i) );
  end
%
%  Compute the Cholesky product.
%
  a(1:n,1:n) = a(1:n,1:n)' * a(1:n,1:n);

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The shifted product U'' * U:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    for j = 1 : n
      fprintf ( 1, '  %14f', a(i,j) );
    end
    fprintf ( 1, '\n' );
  end

  return
end


text:
function [data,interface,options,idata] = loadData(fname,type,varargin)
% import PoleFigure, EBSD, and ODF data
%
% Description
% *loadData* is a low level method for importing EBSD, PoleFigure, ODF and
% Tensor data from external files. It autodetects the format of the file.
% As parameters the method requires a filename and the crystal and specimen
% @symmetry.
%
% Input
%  fname     - filename
%  cs, ss    - crystal, specimen @symmetry (optional)
%
% Options
%  interface  - specific interface to be used
%  comment    - comment to be associated with the data
%
% Output
%  data - @EBSD, @PoleFigure, @SO3Fun, @tensor, @vector3d, @orientation
%
% See also
% ImportEBSDData EBSD/calcODF ebsd_demo loadEBSD_generic

% process input arguments

% get file names
if nargin < 1
  [fname,PathName] = uigetfile( '*.*',...
    'Select Data files');
  fname = [PathName,fname];
end

% read in directory if needed
if ischar(fname)
  if exist(fname,'dir')
    pname = fname;
  else
    pname = fileparts(fname);
  end
  files = dir(fname);
  files = files(~[files.isdir]);
  assert(~isempty(files),'No file found!');
  if ~isempty(pname) && pname(end)~=filesep, pname = [pname,filesep];end
  fname = strcat(pname,{files.name});
end

% get crystal directions
if ~isempty(varargin) && checkClass(varargin{1},'Miller')
  h = vec2cell(varargin{1});
  varargin = varargin(2:end);
end

% get crystal and specimen symmetry
posCS = find(cellfun(@(x) checkClass(x,'symmetry'),varargin),2);

sym = {};
if ~isempty(posCS)
  cs = varargin{posCS(1)};
  if strcmpi(type,'ODF')
    sym = {'cs',cs};
    clear cs;
  end

end

if numel(posCS) ==2
  ss = varargin{posCS(2)};
  if strcmpi(type,'ODF')
    sym = [sym,'ss',{ss}];
    clear ss;
  end
end

varargin(posCS) = [];

% ---------- determine interface ---------------------------

if ~check_option(varargin,'interface')
  [interface,options] = check_interfaces(fname{1},type,varargin{:});
else
  interface = get_option(varargin,'interface');
  options = delete_option(varargin,'interface',1);
end

if isempty(interface), return; end

% --------------- import data ---------------------------------

% determine superposition - PoleFigure only
c = get_option(options,'superposition');
if isa(c,'double'), c = {c};end

if numel(fname) <= 3
  InfoLevel = 1;
else
  InfoLevel = 0;
  hw = waitbar(0,'Loading data files.');
end

for k = 1:numel(fname)
  if exist('hw','var')
    [~,fn,ext] = fileparts(fname{k});
    waitbar(k/numel(fname),hw,['Loading data file ',[fn ext]]);
  end

  data{k} = feval(['load' type '_',char(interface)],...
    fname{k},options{:},sym{:},'InfoLevel',InfoLevel);
end
if exist('hw','var'), close(hw);end

idata = cellfun('prodofsize',data);

% ------------- apply options ----------------------------------

if length(data) == 1 && iscell(data{1}), data = data{1}; end

% set file name
for i = 1:numel(data)
  data{i} = data{i}.setOption('file_name',ls(fname{min(i,length(fname))}));
end

if strcmpi(type,'EBSD') && check_option(varargin,'3d')
  Z = get_option(varargin,'3d',1:numel(data),'double');
  for k=1:numel(data)
    data{k}.z = repmat(Z(k),length(data{k}),1);
  end
  data = [data{:}];
  data.unitCell = calcUnitCell([data.x(:),data.y(:),data.z(:)],varargin{:});
end

% set crystal and specimen symmetry, specimen direction
if ~any(strcmpi(type,{'tensor','vector3d','ODF'}))
  if iscell(data)
    data = cellfun(@(d,f) setOption(d,'file_name',strtrim(ls(f))),data,fname,'UniformOutput',false);
    data = [data{:}];
  end

  if exist('cs','var'), data.CS = cs;end
  if exist('ss','var') && ~isa(data,'EBSD'), data.SS = ss;end % TODO
  if exist('h','var'),  data = set(data,'h',h);end
  if ~isempty_cell(c),  data = set(data,'c',c);end
elseif ~any(strcmpi(type,{'ODF'}))
  data = [data{:}];

  if exist('cs','var')
    if iscell(data)
      for i = 1:numel(data), data{i}.CS = cs; end
    else
      data.CS = cs;
    end
    if exist('ss','var')
      if iscell(data)
        for i = 1:numel(data), data{i}.SS = ss; end
      else
        data.SS = ss;
      end
    end
  end

end

% rotate data
if check_option(varargin,'rotate')
  data = rotate(data,axis2quat(zvector,get_option(varargin,'rotate')));
end

% --------------------------------------------------------------
function v = checkClass(var,className)

if iscell(var) && ~isempty(var)
  v = any(cellfun('isclass',var,className));
else
  v = isa(var,className);
end


text:
function [genes,rxnGeneMat] = getGenesFromGrRules(grRules, originalGenes)
%getGenesFromGrRules  Extract gene list and rxnGeneMat from grRules array.
%
% USAGE:
%
%   [genes,rxnGeneMat] = getGenesFromGrRules(grRules, originalGenes);
%
% INPUTS:
%
%   grRules     A cell array of model grRules, from which a list of genes
%               are to be extracted.
%               NOTE: Boolean operators can be text ("and", "or") or
%                     symbolic ("&", "|"), but there must be a space
%                     between operators and gene names/IDs.
%   originalGenes     The original gene list from the model as reference
%
% OUTPUTS:
%
%   genes       A unique list of all gene IDs that appear in grRules.
%
%   rxnGeneMat  (Optional) A binary matrix indicating which genes
%               participate in each reaction, where rows correspond to
%               reactions (entries in grRules) and columns correspond to
%               genes.
%

% handle input arguments
if nargin < 2
    originalGenes = [];
end

% check if the grRules use written or symbolic boolean operators
if any(contains(grRules,{'&','|'}))
    % fix some potential missing spaces between parentheses and &/|
    grRules = regexprep(grRules,'\)&',') &');   % ")&"  ->  ") &"
    grRules = regexprep(grRules,'&\(','& (');   % "&("  ->  "& ("
    grRules = regexprep(grRules,'\)\|',') |');  % ")|"  ->  ") |"
    grRules = regexprep(grRules,'\|\(','| (');  % "|("  ->  "| ("
else
    % fix some potential missing spaces between parentheses and AND/OR
    grRules = regexprep(grRules,'\)and',') and');  % ")and" ->  ") and"
    grRules = regexprep(grRules,'and\(','and (');  % "and(" ->  "and ("
    grRules = regexprep(grRules,'\)or',') or');    % ")or"  ->  ") or"
    grRules = regexprep(grRules,'or\(','or (');    % "or("  ->  "or ("

    % convert "and" to "&" and "or" to "|" (easier to work with symbols)
    grRules = regexprep(grRules, ' or ', ' | ');
    grRules = regexprep(grRules, ' and ', ' & ');
end

% extract list of genes from each reaction
rxnGenes = cellfun(@(r) unique(regexp(r,'[^&|\(\) ]+','match')),grRules,'UniformOutput',false);

% construct new gene list
nonEmpty = ~cellfun(@isempty,rxnGenes);
genes = unique([rxnGenes{nonEmpty}]');

if ~isempty(originalGenes)
    if ~isequal(sort(originalGenes), sort(genes))
        error('The grRules and original gene list are inconsistent!');
    else
        genes = originalGenes;
    end
end

% construct new rxnGeneMat (if requested)
if nargout > 1
    rxnGeneCell = cellfun(@(rg) ismember(genes,rg),rxnGenes,'UniformOutput',false);
    rxnGeneMat = sparse(double(horzcat(rxnGeneCell{:})'));
end