text:
classdef prtClassMaryLogDisc < prtClass

    properties (SetAccess=private)
        name = 'bleh'  % Logistic Discriminant

        nameAbbreviation = 'meh'  % LogDisc

        isNativeMary = true;  % True
    end

    properties (SetAccess = protected)
        % w
        %   w is a DataSet.nDimensions + 1 x 1 vector of projection weights
        %   learned during LogDisc.train(DataSet)

        wMat = [];  % Regression weights

        % nIterations
        %   Number of iterations used in training.  This is set to a number
        %   between 1 and maxIter during training.

        nIterations = nan;  % The number of iterations used in training
        wChangeTolerance = 1e-2;
        converged = false;
    end

    properties
        % maxIter
        %   Maximum number of iterations to allow before exiting without
        %   convergence.

        maxIter = 500;  % Maxmimuum number of iterations
    end

    methods

        function self = prtClassMaryLogDisc(varargin)
            self = prtUtilAssignStringValuePairs(self,varargin{:});
        end

        function self = set.maxIter(self,val)
            if ~prtUtilIsPositiveScalarInteger(val)
                error('prt:prtClassLogisticDiscriminant:maxIter','maxIter must be a positive scalar integer');
            end
            self.maxIter = val;
        end
    end

    methods (Access=protected, Hidden = true)

        function self = trainAction(self,dataSet)
            %self = trainAction(self,dataSet)

            x = dataSet.getObservations;
            x = cat(2,ones(size(x,1),1),x); %DC component
            y = dataSet.getTargetsAsBinaryMatrix;

            nClasses = dataSet.nClasses;
            d = size(x,2);

            %random initialization; last set of weights is set to 0
            numWeights = d*(nClasses-1);
            weightMatrix = randn(nClasses-1,d);
            weightMatrix = cat(1,weightMatrix,zeros(1,size(weightMatrix,2)));
            weightMatrixOld = weightMatrix;

            %Can calculate B matrix outside loop; makes life fast
            xx = x'*x;
            B = kron(-1/2*(eye(nClasses-1)-ones(nClasses-1)/dataSet.nClasses),xx);

            self.converged = false;
            Binv = B^-1;
            for j = 1:self.maxIter
                psi = (weightMatrix*x')';
                py = bsxfun(@rdivide,exp(psi),sum(exp(psi),2));

                %label error
                yError = y-py;

                %Can we speed this up?  KRON is needlessly slow
                g = 0;
                for i = 1:size(yError,1)
                    g = g + kron(yError(i,:),x(i,:));
                end
                g = g(:);

                % Update weightMatrix in direction of gradient
                weightMatrix(1:end-1,:) = weightMatrix(1:end-1,:) - reshape((Binv*g(1:numWeights)),d,nClasses-1)';

                if norm(weightMatrix(:)-weightMatrixOld(:)) < self.wChangeTolerance
                    self.converged = true;
                    break;
                end
                weightMatrixOld = weightMatrix;

            end
            self.wMat = weightMatrix;
        end

        function ClassifierResults = runAction(self,dataSet)
            %ClassifierResults = runAction(self,DataSet)

            x = dataSet.getObservations;
            x = cat(2,ones(size(x,1),1),x);

            psi = (self.wMat*x')';
            y = bsxfun(@rdivide,exp(psi),sum(exp(psi),2));

            ClassifierResults = dataSet.setObservations(y);
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


%==============================================================================
% This code is part of the Matlab-based toolbox
% LagLDDDM - A Lagrangian Gauss--Newton--Krylov Solver for Mass- and
%                        Intensity-Preserving Diffeomorphic Image Registration
%
% For details and license info see
% - https://github.com/C4IR/FAIR.m/tree/master/add-ons/LagLDDMM
%==============================================================================
%
% Examples for LagLDDMM toolbox.
%
% ELDDMM_2Ddisc2C_mfCurvatureST.m
% ELDDMM_2Ddisc2C_mfDiffusionCC.m
% ELDDMM_2Ddisc2C_mfDiffusionST.m
% ELDDMM_2Ddisc2C_mbCurvatureST.m
% ELDDMM_2Ddisc2C_mbDiffusionCC.m
% ELDDMM_2Ddisc2C_mbDiffusionST.m
% ELDDMM_2Dhands_mfDiffusionCC.m
% ELDDMM_Precond_diffusionCC.m
% ELDDMM_Precond_diffusionST.m
% EMPLDDMM_2DGaussian_mfDiffusionCC.m
% EMPLDDMM_2DGaussian_mfDiffusionST.m
% EMPLDDMM_3Dmouse_mfDiffusionCC.m
% EMPLDDMM_3Dmouse_mfDiffusionST.m
% Ex_HyperElasticNIREP.m
% Ex_LDDMMNIREP.m
%
% ==================================================================================

function debit = contents
if nargout == 0, help(mfilename); return; end;

debit = {
          'ELDDMM_2Ddisc2C_mfCurvatureST.m'
          'ELDDMM_2Ddisc2C_mfDiffusionCC.m'
          'ELDDMM_2Ddisc2C_mfDiffusionST.m'
		  'ELDDMM_2Ddisc2C_mbCurvatureST.m'
          'ELDDMM_2Ddisc2C_mbDiffusionCC.m'
          'ELDDMM_2Ddisc2C_mbDiffusionST.m'
		  'ELDDMM_2Dhands_mfDiffusionCC.m'
          'ELDDMM_Precond_diffusionCC.m'
          'ELDDMM_Precond_diffusionST.m'
          'EMPLDDMM_2DGaussian_mfDiffusionCC.m'
          'EMPLDDMM_2DGaussian_mfDiffusionST.m'
          'EMPLDDMM_3Dmouse_mfDiffusionCC.m'
          'EMPLDDMM_3Dmouse_mfDiffusionST.m'
          'contents.m'
        };


text:
function hline(y, varargin)

% HLINE plot a horizontal line in the current graph

abc = axis;
y = [y y];
x = abc([1 2]);
if length(varargin)==1
  varargin = {'color', varargin{1}};
end
h = line(x, y);
set(h, varargin{:});

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% NTF | bcuNCP | Non-negative CP Decomposition by block-coordinate update (Xu and Yin, 2012)
% process_video('NTF', 'bcuNCP', 'dataset/demo.avi', 'output/bcuNCP.avi');

% Compute a simple NTF model of 10 components
A = double(T);
R = 10; % tensor rank
opts.maxit = 1000; % max number of iterations
opts.tol = 1e-4; % stopping tolerance
M = ncp(T,R,opts);
L = double(full(M));
S = (A - L);

author: andrewssobral

repo: lrslibrary


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

global planC stateS;
indexS = planC{end};

[origin, spacing, center] = getScanOriginSpacing(planC{indexS.scan}(stateS.scanSet));

scanData = planC{indexS.scan}(stateS.scanSet).scanArray;
%scanData = planC{indexS.dose}(stateS.doseSet).doseArray;
% spacing = [1 1 3];

switch upper(command)

    case 'MIP'
        colormap = [];
%         scanData=GPReduce2(scanData,2,0);
        sampleRatio = 2;
        scanData = scanData(1:sampleRatio:end, 1:sampleRatio:end, 1:sampleRatio:end);
        Render3D_MIP(uint16(permute(scanData,[3,1,2])),colormap', 1, [spacing(3),spacing(1),spacing(2)]*10, 1);

    case 'MIP2'
        colormap = [];
        cData3M = scanDoseFusion();
        cData3M = permute(cData3M, [4,1,2,3]);
        [m,n,o,p] = size(cData3M);
        cData3M = reshape(cData3M, [n,o,p*m]);
        alphaRatio = 30;
        Render3D_MIP(uint16(cData3M),colormap', 2, spacing*10, alphaRatio);

    case 'VR1'
%         colormap = [  0.0     0.0   0.0   0.0   0.0;  ...
%                        200.0   0     0     0.8   0.02;  ...
%                        400.0   0.0   1.0   0.0   .8;  ...
%                        600.0   1.0   1.0   0.0   .6; ...
%                        800.0   1.0   0.0   0.0   .7; ...
%                       1000.0   0.0   1.0   1.0   1; ];

        [strVol, colormap] = generateStructVolume();
        sagOn  = 0;
        axOn   = 1;
        corOn  = 1;
        skinOn = 0;
        boneOn = 1;
        sagPos = 130;
        axPos  = 30;
        corPos = 100;
        tRange1 = 2500;
        tRange2 = 1500;
        tRange3 = 1500;
        skinAlpha = 0.3;
        boneAlpha = 0.8;
        Render3D_S(uint16(strVol(:,:,20:end)),colormap', 2, spacing, ...
            sagOn, axOn, corOn, skinOn, boneOn, sagPos, axPos, corPos, tRange1, tRange2, tRange3, skinAlpha, boneAlpha);

    case 'VR2'
        colormap = [   0    1    1   1    .0  ...
                      255   1    1   1    .0  ...
                      800   .9   .9  .9   .1  ...
                      1000  .9   .9  .3   .1  ...
                      1010  1   1  1   1  ...
                      1020  1   1  1   1  ...
                      1030  1   1  1   1  ...
                      1040  1   1  1   1  ...
                      1050  .6   .0  .6   .7  ...
                      1060  .6   .0  .6   .7  ...
                      1070  .6   .0  .6   .7  ...
                      1080  .0   .9  .0   1  ...
                      1090  .0   .6  .0   1  ...
                      1100  .0   .3  .0   1  ...
                      1150  .0   .8  .0   .1  ...
                      1200  .0   .4  .0   .9  ...
                      1500  .0   .6  .0   .9  ...
                      2000  .0   .8  .8   .9  ...
                      4000  1    1   1    1.0 ];

%         scanData = smooth3(scanData, 'gaussian');
        Render3D_V(uint16(scanData(:,:,40:end)),colormap', 2, spacing);

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [growthOnMinimalMedium,essentialExchanges,modelPruned] = predictGrowthRequirements(model)
% Predicts growth requirements for a genome-scale reconstruction.
% The list of predicted essential exchanges as well as a model
% constrained with these exchanges are returned.
%
% INPUT
% model                  COBRA model structure
% biomassReaction        String listing the biomass reaction
%
% OUTPUT
% growthOnMinimalMedium  Bool if growth on defined medium yes or no
% essentialExchanges     Exchanges that need to be open to enable growth
% constrainedModel       Model constrained with predicted essential
%                        exchanges
%
% Almut Heinken, November 2020

tol=1e-6;

% set "unlimited" constraints
model = changeRxnBounds(model, model.rxns(strncmp('EX_', model.rxns, 3)), -1000, 'l');
model = changeRxnBounds(model, model.rxns(strncmp('EX_', model.rxns, 3)), 1000, 'u');

% remove positive values on lower bounds
model.lb(find(model.lb>0))=0;
%
% % remove compounds that should not be essential for any organism
% unlikelyRequirements={'EX_alaasp(e)';'EX_alagln(e)';'EX_alaglu(e)';'EX_alagly(e)';'EX_alahis(e)';'EX_alaleu(e)';'EX_alathr(e)';'EX_glygln(e)';'EX_glyglu(e)';'EX_glyleu(e)';'EX_glymet(e)';'EX_glyphe(e)';'EX_glypro(e)';'EX_glytyr(e)';'EX_metala(e)';'EX_glyasn(e)';'EX_glyasp(e)';'EX_dpcoa(e)';'EX_coa(e)';'EX_cmp(e)';'EX_dtmp(e)';'EX_ump(e)';'EX_gmp(e)';'EX_dgtp(e)';'EX_nadp(e)';'EX_nad(e)';'EX_thmmp(e)';'EX_malt(e)';'EX_stys(e)';'EX_raffin(e)';'EX_melib(e)';'EX_lcts(e)';'EX_gtp(e)';'EX_datp(e)';'EX_amp(e)';'EX_pppi(e)';'EX_r5p(e)';'EX_g1p(e)';'EX_glu_D(e)';'EX_agm(e)';'EX_glycys(e)';'EX_1hibup_S(e)';'EX_1hibupglu_S(e)';'EX_1hmdgluc(e)';'EX_1ohmdz(e)';'EX_2hatvacid(e)';'EX_2hatvacidgluc(e)';'EX_2hatvlac(e)';'EX_2hatvlacgluc(e)';'EX_2hibup_S(e)';'EX_2hibupglu_S(e)';'EX_2oh_cbz_glc(e)';'EX_2oh_cbz(e)';'EX_2oh_mtz_glc(e)';'EX_2oh_mtz(e)';'EX_34dhphe(e)';'EX_3hibup_S(e)';'EX_3hibupglu_S(e)';'EX_3meacmp(e)';'EX_3oh_cbz_glc(e)';'EX_3oh_cbz(e)';'EX_3oh_dlor_glc(e)';'EX_3oh_dlor(e)';'EX_3oh_mdea_glc(e)';'EX_3oh_mdea(e)';'EX_3oh_mxn_glc(e)';'EX_3oh_mxn(e)';'EX_4dh_tpno_1glc(e)';'EX_4dh_tpno(e)';'EX_4hmdgluc(e)';'EX_4hphac(e)';'EX_4oh_dcf_glc(e)';'EX_4oh_dcf(e)';'EX_4oh_kp_glc(e)';'EX_4oh_kp(e)';'EX_4oh_levole_glc(e)';'EX_4oh_levole(e)';'EX_4oh_meth_glc(e)';'EX_4oh_meth(e)';'EX_4oh_propl_glc(e)';'EX_4oh_propl(e)';'EX_4oh_trz_glc(e)';'EX_4oh_trz(e)';'EX_4oh_vcz_glc(e)';'EX_4oh_vcz(e)';'EX_4ohmdz(e)';'EX_5asa(e)';'EX_5fura(e)';'EX_5oh_sulfp_glc(e)';'EX_5oh_sulfp(e)';'EX_5ohfvs(e)';'EX_5ohfvsglu(e)';'EX_6bhglz(e)';'EX_6bhglzglc(e)';'EX_6ohfvs(e)';'EX_6ohfvsglu(e)';'EX_7a_czp(e)';'EX_7bhglz(e)';'EX_7bhglzglc(e)';'EX_7oh_efv_glc(e)';'EX_7oh_efv(e)';'EX_814dioh_efv_glc(e)';'EX_814dioh_efv(e)';'EX_8oh_efv_glc(e)';'EX_8oh_efv(e)';'EX_abz_ala_b(e)';'EX_ac_amn_b_gly_glc(e)';'EX_ac_amn_b_gly(e)';'EX_ac5asa(e)';'EX_acisnzd(e)';'EX_acmp_glc(e)';'EX_acmp(e)';'EX_acmpglu(e)';'EX_alpz_4oh_glc(e)';'EX_alpz_4oh(e)';'EX_alpz_aoh_glc(e)';'EX_alpz_aoh(e)';'EX_am14_glc(e)';'EX_am14(e)';'EX_am1ccs(e)';'EX_am1cglc(e)';'EX_am5_glc(e)';'EX_am5(e)';'EX_am6_glc(e)';'EX_am6(e)';'EX_amio_c_glc(e)';'EX_amio_c(e)';'EX_amio_glc(e)';'EX_amio(e)';'EX_amn_b_gly_glc(e)';'EX_amn_b_gly(e)';'EX_amntd_m6_glc(e)';'EX_amntd_m6(e)';'EX_anzp(e)';'EX_atvacid(e)';'EX_atvacylgluc(e)';'EX_atvethgluc(e)';'EX_atvlac(e)';'EX_atvlacgluc(e)';'EX_bhpm_glc(e)';'EX_bhpm(e)';'EX_bilr_355(e)';'EX_bilr_M10(e)';'EX_bilr_M12(e)';'EX_bilr_M14(e)';'EX_bilr_M15(e)';'EX_bilr_M16(e)';'EX_bilr_M7(e)';'EX_brv(e)';'EX_bsn_glc(e)';'EX_bsn(e)';'EX_bvu(e)';'EX_bz_glc(e)';'EX_bz(e)';'EX_bzd(e)';'EX_C02528(e)';'EX_caribup_s(e)';'EX_caribupglu_S(e)';'EX_cbz_glc(e)';'EX_cbz(e)';'EX_cd6168_glc(e)';'EX_cd6168(e)';'EX_chlphncl_glc(e)';'EX_chlphncl(e)';'EX_cholate(e)';'EX_clcxb_c_glc(e)';'EX_clcxb_c(e)';'EX_clcxb_glc(e)';'EX_clcxb(e)';'EX_clobi_c(e)';'EX_clobi_glc(e)';'EX_crvsm1(e)';'EX_crvsm23(e)';'EX_cvm1gluc(e)';'EX_cvm23gluc(e)';'EX_czp(e)';'EX_daa_glc(e)';'EX_daa(e)';'EX_dcf_glc(e)';'EX_dcf(e)';'EX_dchac(e)';'EX_ddea_glc(e)';'EX_ddea(e)';'EX_des_astzl_glc(e)';'EX_des_astzl(e)';'EX_dfdcytd(e)';'EX_dfduri(e)';'EX_dgchol(e)';'EX_dh5fura(e)';'EX_digitoxin(e)';'EX_digoxin_glc(e)';'EX_digoxin(e)';'EX_dihydro_digitoxin(e)';'EX_dihydro_digoxin(e)';'EX_dlb_glc(e)';'EX_dlb(e)';'EX_dnpz_5des(e)';'EX_dnpz_6des(e)';'EX_dnpz_m11(e)';'EX_dnpz_m12(e)';'EX_dnpz_m13(e)';'EX_dnpz_m14(e)';'EX_dnpz_m9(e)';'EX_doh_etr_glc(e)';'EX_doh_etr(e)';'EX_doh_vcz_glc(e)';'EX_doh_vcz(e)';'EX_dopa(e)';'EX_dxo_glc(e)';'EX_dxo(e)';'EX_efv_glc(e)';'EX_efv(e)';'EX_eltr_glc(e)';'EX_eltr_m3(e)';'EX_eltr_m4(e)';'EX_eltr(e)';'EX_estradiol(e)';'EX_estradiolglc(e)';'EX_estriol(e)';'EX_estriolglc(e)';'EX_estrone(e)';'EX_estroneglc(e)';'EX_eztmb_glc(e)';'EX_eztmb(e)';'EX_fcsn(e)';'EX_fvs(e)';'EX_fvsgluc(e)';'EX_fvstet(e)';'EX_fvstetglu(e)';'EX_gchola(e)';'EX_glc3meacp(e)';'EX_gltmn_glc(e)';'EX_gltmn(e)';'EX_gmfl_glc(e)';'EX_gmfl_mI_glc(e)';'EX_gmfl_mI(e)';'EX_gmfl_mII_glc(e)';'EX_gmfl_mII(e)';'EX_gmfl_mIII_glc(e)';'EX_gmfl_mIII(e)';'EX_gmfl(e)';'EX_gtacmp(e)';'EX_HC02191(e)';'EX_hst_3_glc(e)';'EX_hst_3_s(e)';'EX_hst_37_diglc(e)';'EX_hst_3glc_7s(e)';'EX_hst_7_glc(e)';'EX_hst_7_s(e)';'EX_hst_7glc_3s(e)';'EX_hst(e)';'EX_ibup_S(e)';'EX_ibupgluc(e)';'EX_imn_glc(e)';'EX_imn(e)';'EX_inv_m1(e)';'EX_inv(e)';'EX_isnzd(e)';'EX_isosorbide_5mn_glc(e)';'EX_isosorbide_5mn(e)';'EX_kprofen_glc(e)';'EX_kprofen(e)';'EX_lactl(e)';'EX_lst4exp(e)';'EX_lstn(e)';'EX_lstn1gluc(e)';'EX_lstnm4(e)';'EX_lstnm7(e)';'EX_mdz_glc(e)';'EX_mdz(e)';'EX_mdzglc(e)';'EX_miso_glc(e)';'EX_miso(e)';'EX_mrphn_3glc(e)';'EX_mrphn_6glc(e)';'EX_mrphn(e)';'EX_mtym(e)';'EX_mtz_glc(e)';'EX_mtz(e)';'EX_N_oh_phtn_glc(e)';'EX_N_oh_phtn(e)';'EX_nchlphncl(e)';'EX_neopront(e)';'EX_nsldp_m5_glc(e)';'EX_nsldp_m5(e)';'EX_nverp_glc(e)';'EX_nverp(e)';'EX_nzp(e)';'EX_odsm_egltmn_glc(e)';'EX_odsm_egltmn(e)';'EX_odsm_gltmn_glc(e)';'EX_odsm_gltmn(e)';'EX_oh_etr_glc(e)';'EX_oh_etr(e)';'EX_oh_pbl_glc(e)';'EX_oh_pbl(e)';'EX_olsa(e)';'EX_pcresol(e)';'EX_phppa_glc(e)';'EX_phppa(e)';'EX_phtn_glc(e)';'EX_prob_glc(e)';'EX_prob(e)';'EX_pront_glc(e)';'EX_pront(e)';'EX_propl_glc(e)';'EX_propl(e)';'EX_prx_mI_glc(e)';'EX_prx_mI(e)';'EX_prx_mII_glc(e)';'EX_prx_mII(e)';'EX_ptvst(e)';'EX_ptvstgluc(e)';'EX_pvs(e)';'EX_pvsgluc(e)';'EX_R_6oh_warf_glc(e)';'EX_R_6oh_warf(e)';'EX_R_7oh_warf_glc(e)';'EX_R_7oh_warf(e)';'EX_R_8oh_warf_glc(e)';'EX_R_8oh_warf(e)';'EX_r406_glc(e)';'EX_r406(e)';'EX_r529_glc(e)';'EX_r529(e)';'EX_r788(e)';'EX_regfnb_glc(e)';'EX_regfnb(e)';'EX_rep_glc(e)';'EX_rep(e)';'EX_rpn_104557_cb_glc(e)';'EX_rpn_104557(e)';'EX_rpn_96990_glc(e)';'EX_rpn_96990(e)';'EX_rpn_oh_glc(e)';'EX_rpn_oh(e)';'EX_rsv(e)';'EX_rsvgluc(e)';'EX_S_4oh_warf_glc(e)';'EX_S_4oh_warf(e)';'EX_S_6oh_warf_glc(e)';'EX_S_6oh_warf(e)';'EX_sanilamide(e)';'EX_sb_611855(e)';'EX_sb_y(e)';'EX_sch_488128(e)';'EX_sch_57871_glc(e)';'EX_sch_57871(e)';'EX_sfnd_1689_glc(e)';'EX_sfnd_1689(e)';'EX_sftz_glc(e)';'EX_sftz(e)';'EX_simvgluc(e)';'EX_smap_glc(e)';'EX_smap(e)';'EX_smvacid(e)';'EX_sn38(e)';'EX_sn38g(e)';'EX_spz_glc(e)';'EX_spz_sfn_glc(e)';'EX_spz_sfn(e)';'EX_spz(e)';'EX_srv(e)';'EX_ssz(e)';'EX_stg_m3(e)';'EX_stg_m4(e)';'EX_stg(e)';'EX_sulfp(e)';'EX_tab(e)';'EX_tat(e)';'EX_taur(e)';'EX_tchola(e)';'EX_tdchola(e)';'EX_tgz_glc(e)';'EX_tgz(e)';'EX_thsacmp(e)';'EX_tlf_a_1a(e)';'EX_tlf_a_1b(e)';'EX_tlf_a_1x(e)';'EX_tlf_a_2(e)';'EX_tlf_a_2a(e)';'EX_tlf_a_3(e)';'EX_tlf_a_4(e)';'EX_tlf_a_m1(e)';'EX_tlf_a_m2(e)';'EX_tlf_a_m3(e)';'EX_tlf_a_m4(e)';'EX_tlf_a_m4a(e)';'EX_tlf_a_m5(e)';'EX_tlf_a_m5a(e)';'EX_tlf_a_m5b(e)';'EX_tlf_a_m6(e)';'EX_tlf_a_m9(e)';'EX_tlf_a(e)';'EX_tlms_glc(e)';'EX_tlms(e)';'EX_tmacmp(e)';'EX_tolcp_ac_glc(e)';'EX_tolcp_ac(e)';'EX_tolcp_am_glc(e)';'EX_tolcp_am(e)';'EX_tolcp_glc(e)';'EX_tolcp(e)';'EX_tpno_1glc_4g(e)';'EX_tpno_4g(e)';'EX_tpno_4glc(e)';'EX_tpnoh(e)';'EX_tsacmgluc(e)';'EX_tym(e)'};
% model=changeRxnBounds(model,unlikelyRequirements,0,'l');

% set objective
biomassReaction=model.rxns{find(strncmp(model.rxns,'bio',3)),1};
model = changeObjective(model, biomassReaction);

[modelMin, modelPruned, essentialExchanges] = generateCompactExchModel(model,tol,biomassReaction,1,1);
modelMin = changeObjective(modelMin, biomassReaction);
FBA=optimizeCbModel(modelMin,'max');
growthOnMinimalMedium=FBA.f;

end

author: opencobra

repo: cobratoolbox


text:
function h = showmesh(node,elem,varargin)
%% SHOWMESH displays a triangular mesh in 2-D.
%
%    showmesh(node,elem) displays a topological 2-dimensional mesh,
%    including planar meshes and surface meshes. The mesh is given by node
%    and elem matrices; see meshdoc for the mesh data structure:
%    node and elem.
%
%    showmesh(node,elem,viewangle) changes the display angle. The
%    deault view angle on planar meshes is view(2) and view(3) for surface
%    meshes.
%
%    showmesh(node,elem,'param','value','param','value'...) allows
%    additional patch param/value pairs to be used when displaying the
%    mesh.  For example, the default transparency parameter for a surface
%    mesh is set to 0.75. You can overwrite this value by using the param
%    pair ('FaceAlpha', value). The value has to be a number between 0 and
%    1. Other parameters include: 'Facecolor', 'Edgecolor' etc. These
%    parameters are mostly used when displaying a surface mesh.
%
%    To display a 3-dimensional mesh, use showmesh3 or showboundary3.
%
%   Example:
%     % A mesh for a L-shaped domain
%     [node,elem] = squaremesh([-1,1,-1,1],0.5);
%     [node,elem] = delmesh(node,elem,'x>0 & y<0');
%     figure;
%     showmesh(node,elem);
%
%     % A mesh for a unit sphere
%     node = [1,0,0; 0,1,0; -1,0,0; 0,-1,0; 0,0,1; 0,0,-1];
%     elem = [6,1,2; 6,2,3; 6,3,4; 6,4,1; 5,1,4; 5,3,4; 5,3,2; 5,2,1];
%     for i = 1:3
%         [node,elem] = uniformrefine(node,elem);
%     end
%     r = sqrt(node(:,1).^2 + node(:,2).^2 + node(:,3).^2);
%     node = node./[r r r];
%     figure;
%     subplot(1,2,1);
%     showmesh(node,elem);
%     subplot(1,2,2);
%     showmesh(node,elem,'Facecolor','y','Facealpha',0.5);
%
%   See also showmesh3, showsolution, showboundary3.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

dim = size(node,2);
nv = size(elem,2);
if (dim==2) && (nv==3) % planar triangulation
    h = trisurf(elem(:,1:3),node(:,1),node(:,2),zeros(size(node,1),1));
    set(h,'facecolor',[0.5 0.9 0.45],'edgecolor','k');
    view(2); axis equal; axis tight; axis off;
end
if (dim==2) && (nv==4) % planar quadrilateration
    h = patch('Faces', elem, 'Vertices', node);
    set(h,'facecolor',[0.5 0.9 0.45],'edgecolor','k');
    view(2); axis equal; axis tight; axis off;
end
if (dim==3)
    if size(elem,2) == 3 % surface meshes
        h = trisurf(elem(:,1:3),node(:,1),node(:,2),node(:,3));
        set(h,'facecolor',[0.5 0.9 0.45],'edgecolor','k','FaceAlpha',0.75);
        view(3); axis equal; axis off; axis tight;
    elseif size(elem,3) == 4
        showmesh3(node,elem,varargin{:});
        return
    end
end
if (nargin>2) && ~isempty(varargin) % set display property
    if isnumeric(varargin{1})
        view(varargin{1});
        if nargin>3
            set(h,varargin{2:end});
        end
    else
        set(h,varargin{1:end});
    end
end

repo: ifem


text:
classdef MMF7 < PROBLEM
% <multi> <real> <multimodal>
% Multi-modal multi-objective test function

%------------------------------- Reference --------------------------------
% C. Yue, B. Qu, and J. Liang, A multi-objective particle swarm optimizer
% using ring topology for solving multimodal multiobjective Problems, IEEE
% Transactions on Evolutionary Computation, 2018, 22(5): 805-817.
%------------------------------- 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
        POS;    % Pareto optimal set for IGDX calculation
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            obj.D = 2;
            obj.lower    = [1,-1];
            obj.upper    = [3,1];
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            PopObj(:,1) = abs(X(:,1)-2);
            PopObj(:,2) = 1 - sqrt(PopObj(:,1)) + (X(:,2)-(0.3.*(PopObj(:,1).^2).*cos(24.*pi.*PopObj(:,1)+4.*pi)+0.6.*PopObj(:,1)).*sin(6.*pi.*PopObj(:,1)+pi)).^2;
        end
        %% Generate Pareto optimal solutions
        function R = GetOptimum(obj,N)
            % Generate points in Pareto optimal set
            obj.POS(:,1) = linspace(1,3,N)';
            obj.POS(:,2) = (0.3*(obj.POS(:,1)-2).^2.*cos(24*pi*abs(obj.POS(:,1)-2)+4*pi)+0.6*abs(obj.POS(:,1)-2)).*sin(6*pi*abs(obj.POS(:,1)-2)+pi);
            % Generate points on Pareto front
            R(:,1) = linspace(0,1,N)';
            R(:,2) = 1 - sqrt(R(:,1));
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R(:,1) = linspace(0,1,100)';
            R(:,2) = 1 - sqrt(R(:,1));
        end
        %% Calculate the metric value
        function score = CalMetric(obj,metName,Population)
            switch metName
                case 'IGDX'
                    score = feval(metName,Population,obj.POS);
                otherwise
                    score = feval(metName,Population,obj.optimum);
            end
        end
        %% Display a population in the objective space
        function DrawObj(obj,Population)
            PopDec = Population.decs;
            temp   = PopDec(:,1) <= 2;
            Draw(Population(temp).objs,'o','MarkerSize',6,'Marker','o','Markerfacecolor',[1 .5 .5],'Markeredgecolor',[1 .2 .2],{'\it f\rm_1','\it f\rm_2',[]});
            Draw(Population(~temp).objs+0.1,'o','MarkerSize',6,'Marker','o','Markerfacecolor',[.5 .5 1],'Markeredgecolor',[.2 .2 1]);
            Draw(obj.PF,'-','LineWidth',1,'Color',[1 .2 .2]);
            Draw(obj.PF+0.1,'-','LineWidth',1,'Color',[.2 .2 1]);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
 function x = dtft_adj_arrayfun(X, omega, Nd, n_shift, usearrayfun)
%function x = dtft_adj(X, omega, Nd, n_shift, usecellfun)
%|
%| Compute adjoint of d-dim DTFT for spectrum X at frequency locations omega
%|
%| in
%|	X	[M L]		dD DTFT values
%|	omega	[M d]		frequency locations (radians)
%|	n_shift [d 1]		use [0:N-1]-n_shift (default [0 ... 0])
%|	usearrayfun			1 to reduce memory use (slower)
%| out
%|	x	[(Nd) L]	signal values
%|
%| Requires enough memory to store M * (*Nd) size matrices. (For testing only.)
%|
%| Copyright 2003-4-13, Jeff Fessler, University of Michigan
%| Revised 2013-3-22, Daniel Weller, University of Michigan

% if no arguments, then run a simple test
if nargin < 2
	help(mfilename)
	Nd = [4 6 5]; Nd = [50 20 10];
	n_shift = [2 1 3];
	n_shift = 0*[2 1 3]; n_shift = [20 10 4];
	% test with uniform frequency locations:
	o1 = 2*pi*[0:(Nd(1)-1)]'/Nd(1);
	o2 = 2*pi*[0:(Nd(2)-1)]'/Nd(2);
	o3 = 2*pi*[0:(Nd(3)-1)]'/Nd(3);
	[o1 o2 o3] = ndgrid(o1, o2, o3);
	X = o1 + o2 - o3; % test spectrum
	om = [o1(:) o2(:) o3(:)];
	td = tic(); xd = dtft_adj(X(:), om, Nd, n_shift); td = toc(td);
	tl = tic(); xl = dtft_adj(X(:), om, Nd, n_shift, 1); tl = toc(tl);
	printm('loop max %% difference = %g (te = %g/%g)', max_percent_diff(xl,xd), tl, td)
	tc = tic(); xc = dtft_adj_arrayfun(X(:), om, Nd, n_shift, 1); tc = toc(tc);
	printm('arrayfun max %% difference = %g (te = %g/%g)', max_percent_diff(xc,xd), tc, td)
	Xp = X .* reshape(exp(-1i * om * n_shift(:)), size(X));
	xf = ifftn(Xp) * prod(Nd);
	printm('ifft max %% difference = %g', max_percent_diff(xf,xd))
return
end

if ~isvar('n_shift') || isempty(n_shift), n_shift = zeros(size(Nd)); end
if ~isvar('usearrayfun') || isempty(usearrayfun), usearrayfun = 0; end

% if length(Nd) == 1
%	nn{1} = [0:(Nd(1)-1)] - n_shift(1);
% elseif length(Nd) == 2
%	nn{1} = [0:(Nd(1)-1)] - n_shift(1);
%	nn{2} = [0:(Nd(2)-1)] - n_shift(2);
%	[nn{1} nn{2}] = ndgrid(nn{1}, nn{2});
% elseif length(Nd) == 3
%	nn{1} = [0:(Nd(1)-1)] - n_shift(1);
%	nn{2} = [0:(Nd(2)-1)] - n_shift(2);
%	nn{3} = [0:(Nd(3)-1)] - n_shift(3);
%	[nn{1} nn{2} nn{3}] = ndgrid(nn{1}, nn{2}, nn{3});
% else
%	'only 1D-3D done'
% end
dd = size(omega,2);
Nd(end+1:dd) = 1;
n_shift(end+1:dd) = 0;

nn = cell(1,dd);
for id=1:dd % fixed: dd
	nn{id} = (0:(Nd(id)-1))-n_shift(id);
end
[nn{:}] = ndgrid(nn{:}); % fixed: dd
nn = cellfun(@(x) col(x),nn,'UniformOutput',false);
nn = cat(2,nn{:}); % [*Nd dd]
omega = omega.'; % [dd M]

if usearrayfun
    x = arrayfun(@(n) exp(1i.*(nn(n,:)*omega)) * X,1:prod(Nd),'UniformOutput',false);
    x = cat(1,x{:});
else
    x = exp(1i.*(nn*omega)) * X;
end
x = reshape(x, [Nd ncol(X)]); % [Nd L]

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function f = ddivWexp(Wxk_r, varargin)
% function f = ddivW(Wxk_r, varargin)
% Wxk_r = reshape(Wxk,1,x*k) -> row vector
% Vxt = varargin{1};        %data
% Hkt = varargin{2};        %H matrix
% Wxk_fix=varargin{3};      %fixed part of the Wxk matrix (e.g. background)
% Hkt_fix = varargin{4};    %fixed part (lines) of the H matrix (e.g. background)

Vxt = varargin{1};      %data
Hkt_tmp = varargin{2};  %H matrix
Wxk_fix = varargin{3};  %fixed part of the Wxk matrix (e.g. background) ->rows
Hkt_fix = varargin{4};  %fixed part (lines) of the H matrix (e.g. background)
peval = varargin{5}; %parameters

if ~isfield(peval, 'w_lambda') peval.w_lambda=0; end

x=size(Vxt,1);
k=length(peval.w_dovec);

Wxk_tmp = exp(reshape(Wxk_r,x,k));

Wxk = zeros(peval.numpix, peval.ncomp+1);
Hkt = zeros(peval.ncomp+1, peval.nt);

Wxk(:,peval.w_dovec)=Wxk_tmp;
Hkt(peval.h_dovec,:)=Hkt_tmp;

Wxk(:,peval.w_fixvec)=Wxk_fix;
Hkt(peval.h_fixvec,:)=Hkt_fix;

%%%%
% sumw = sum(Wxk,1);
% Wxk = Wxk./repmat(sumw,size(Wxk,1),1);

% fxt = Vxt.*log(Vxt./(Wxk*Hkt))-Vxt+Wxk*Hkt; %d-divergence
% f = sum(fxt(:));
% peval.w_lambda*Wxk.^2
f = ddivergence(Vxt,Wxk*Hkt);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


%  Script file: test_preallocate.m
%
%  Purpose:
%    This program tests the creation of cell arrays with and
%    without preallocation.
%
%  Record of revisions:
%      Date       Programmer          Description of change
%      ====       ==========          =====================
%    03/04/07    S. J. Chapman        Original code
%
% Define variables:
%   a            -- Cell array
%   maxvals      -- Maximum values in cell array

% Create array without preallocation
clear all
maxvals = 50000;
tic
for ii = 1:maxvals
   a{ii} = ['Element ' int2str(ii)];
end
disp( ['Elapsed time without preallocation = ' num2str(toc)] );

% Create array with preallocation
clear all
maxvals = 50000;
tic
a = cell(1,maxvals);
for ii = 1:maxvals
   a{ii} = ['Element ' int2str(ii)];
end
disp( ['Elapsed time with preallocation    = ' num2str(toc)] );


text:
function [] = get_quiver_data

ax = gca;
h = get(ax, 'children');
hq = findobj(h, 'flat', 'type', 'hggroup');

q = get(hq);

c = get(hq, 'children');

get(c(1) )
xb = get(c(1), 'XData');
yb = get(c(1), 'YData');
zb = get(c(1), 'ZData');

ax1 = newax;

plotmd(ax1, [xb; yb; zb] )

hold(ax1, 'on')

get(c(1) )
xh = get(c(2), 'XData');
yh = get(c(2), 'YData');
zh = get(c(2), 'ZData');

plotmd(ax1, [xh; yh; zh] )
axis equal

author: Sable

repo: mcbench-benchmarks


text:
function K = k00(x, xp, hyp, ubar, ubarp, dt, i)

K = k00.ku0u0(x, xp, hyp, ubar, ubarp, dt, i);

end

author: maziarraissi

repo: HPM

save_path: github-repos/MATLAB/maziarraissi-HPM


text:
function toms462_test02 ( )

%*****************************************************************************80
%
%% TOMS462_TEST02 tests BIVNOR.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 May 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TOMS462_TEST02\n' );
  fprintf ( 1, '  Compare BIVNOR with some tabulated data.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X          Y          ' );
  fprintf ( 1, 'R           P                         P' );
  fprintf ( 1, '                      DIFF\n' );
  fprintf ( 1, '                                ' );
  fprintf ( 1, '       (Tabulated)               (BIVNOR)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, y, r, fxy1 ] = bivariate_normal_cdf_values ( n_data );

    if ( n_data == 0 )
      break
    end
%
%  BIVNOR computes the "tail" of the probability, and we want the
%  initial part so take a short cut and just pass in the negatives of X and Y.
%
    fxy2 = bivnor ( - x, - y, r );

    fprintf ( 1, '  %8.4f  %8.4f  %8.4f  %24.16g  %24.16g  %10.4g\n', ...
       x, y, r, fxy1, fxy2, abs ( fxy1 - fxy2 ) );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
mex ind2subv.c
mex subv2ind.c
mex normalise.c
mex -c mexutil.c
if ~isunix
  mex repmatC.c mexutil.obj
else
  mex repmatC.c mexutil.o
end


text:
%Run one participant: normalize full brain tractography set (not just mori).
%Two-stage procedure: first, generate a mori set and initialize the solution. Second, throw in the rest of the fibers and normalize them, too.

%Assumes each dti folder already contains safe brain mask, bvals/bvecs averages & contrack parameter files.
%See dti_Longitude_NormalizeMoriFibersDensity.m for an example how to set up a batch process to generate those.
%In face, use dtiRawAverage instead of dtiAverageRawAligned function (the fist one is a more general version).

%ER 07/2009
warning('off'); clear;
addpath(genpath('~/vistasoft'));
addpath(genpath('/usr/local/matlab/toolbox/mri/spm5_r2008/'));

project_folder='/biac3/wandell4/data/reading_longitude/dti_y1234/';
load('/biac3/wandell4/users/elenary/longitudinal/subjectCodesAll4Years');
subjectID=subjectCodes;
numProcessorsAvail=16; %Courtesy max number of processors to be used on Azure without annoying everyone else.
allfibers_prefix='all';
moriFile=fullfile('dti06trilinrt', 'fibers',  [allfibers_prefix 'ConnectingGM_MoriGroups.mat']);
wholebrainGM=fullfile('dti06trilinrt', 'fibers',  [allfibers_prefix 'ConnectingGM.mat']);
wholebrainFibers=fullfile('dti06trilinrt', 'fibers', [allfibers_prefix '.mat']); %Usually we do not save all.mat since blue
%matter only work with fibers that connect gm

compute_allConnectingGM=false; %We already computed allConnectingGM before; only recompute if does not exist
saveAllFibers=false;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for s=1:size(subjectID, 2)

    cd([project_folder subjectID{s}]);
    if exist(moriFile, 'file')
    fprintf('Subject %s already has morifile created \n', num2str(s));
    continue;
    end

    fprintf(1, 'Working on %s\n', subjectID{s});
    dt6File=fullfile(pwd, 'dti06trilinrt', 'dt6.mat');
    dt = dtiLoadDt6(dt6File);

    if  compute_allConnectingGM || ~exist(wholebrainGM, 'file')
    %Perform full brain tractography and save results
        faThresh = 0.30;
        opts.stepSizeMm = 1;
        opts.faThresh = 0.15;
        opts.lengthThreshMm = [50 250];
        opts.angleThresh = 50;
        opts.wPuncture = 0.2;
        opts.whichAlgorithm = 1;
        opts.whichInterp = 1;
        opts.seedVoxelOffsets = [0.25 0.75];
        opts.offsetJitter = 0.1;
        fa = dtiComputeFA(dt.dt6);
        fa(fa>1) = 1; fa(fa<0) = 0;
        roiAll = dtiNewRoi(allfibers_prefix);
        mask = dtiCleanImageMask(fa>=faThresh);
        [x,y,z] = ind2sub(size(mask), find(mask));
        clear mask fa;
        roiAll.coords = mrAnatXformCoords(dt.xformToAcpc, [x,y,z]);
        clear x y z;
        fg = dtiFiberTrack(dt.dt6, roiAll.coords, dt.mmPerVoxel, dt.xformToAcpc, allfibers_prefix, opts);
        clear roiAll
    if saveAllFibers
        dtiWriteFiberGroup(fg, wholebrainFibers);
    end
    fprintf('Get fibers connecting GM\n');
    fg=dtiGetFibersConnectingGM(fg, dt);
    dtiWriteFiberGroup(fg, wholebrainGM); %Need to save it once b/c it is fed to dtiFindMoriTracts; overwritten later by a sequence of Mori and Non-Mori (unclassified) fibers.
    end

    fprintf('Find Mori Tracts\n');

    [fgMoriWithCST, fg_unclassifiedMoriWithCST]=dtiFindMoriTracts(dt6File, moriFile, wholebrainGM);%Note: here we went with default values which included useRoiBasedApproach=[.89 1]. This is a very conservative threshold for minDist (fiber-to-roi). But since we are simply using this moriClassification to generate an initial solution for blue Matter, we'll leave it at that.
    mergedFG = dtiMergeFiberGroups(fgMoriWithCST, fg_unclassifiedMoriWithCST, allfibers_prefix);
    %Save "Mori-plus-the-rest" file (overwrites/reorders wholebrainGM)
    dtiWriteFiberGroup(mergedFG, wholebrainGM);
    %Save Mori only file
    dtiWriteFiberGroup(fgMoriWithCST, moriFile);
    %Save both of them as PDB
    dtiWriteFibersPdb(mergedFG, dt.xformToAcpc, [wholebrainGM(1:(end-4)) '.pdb']);
    dtiWriteFibersPdb(fgMoriWithCST, dt.xformToAcpc, [moriFile(1:(end-4)) '.pdb']);

    %Create initial solution
    [status, result]=system(['ps -U elenary |grep ''trueSA'' |wc|awk ''{print $1}''']); %result=2 if no trueSA by me are running
    if str2num(result)<=(numProcessorsAvail+2)
    display(['starting a trueSA for' subjectID{s}]);
    SAlogfile=['/azure/scr1/bluematterscratch' filesep subjectID{s} 'dti06trilinrt_' allfibers_prefix 'ConnectingGM_MoriGroups_DN' filesep 'runSAmpi.log'];
    system(['/biac3/wandell4/users/elenary/density_normalization/scripts/runSAmpi.sh ' project_folder filesep subjectID{s} filesep 'dti06trilinrt  ' [allfibers_prefix 'ConnectingGM_MoriGroups.pdb'] ' ' [allfibers_prefix 'ConnectingGM_MoriGroups_DN.pdb'] ' &>' SAlogfile  '&']);
    else continue
    end
end
return
%Since we could not launch all trueSAs, launch'em here
a=1;
while a
try
    a=(size(findstr('allConnectingGM_MoriGroups_DN.ind', ls([project_folder filesep '*' filesep 'dti06trilinrt' filesep 'fibers' filesep allfibers_prefix  'ConnectingGM_MoriGroups_DN.ind'])), 2)<108);
catch
    %in no *.ind yet
    a=1;
end
for s=1:size(subjectID, 2)
    s
     cd([project_folder subjectID{s}]);
        [status, result]=system('ps -U elenary -u elenary |grep ''trueSA'' |wc|awk ''{print $1}''');
    if (~exist([moriFile(1:(end-4)) '_DN.TMP'], 'file')) && str2num(result)<=(numProcessorsAvail+2) && (~exist([moriFile(1:(end-4)) '_DN.ind'], 'file')) && (exist(moriFile, 'file'))
        display(['starting a trueSA for ' subjectID{s}]);
        SAlogfile=['/azure/scr1/bluematterscratch' filesep subjectID{s} 'dti06trilinrt_' allfibers_prefix 'ConnectingGM_MoriGroups_DN' filesep 'runSAmpi.log'];
    mkdir(fileparts(SAlogfile));
        system(['/biac3/wandell4/users/elenary/density_normalization/scripts/runSAmpi.sh ' project_folder filesep subjectID{s} filesep 'dti06trilinrt  ' [allfibers_prefix 'ConnectingGM_MoriGroups.pdb'] ' ' [allfibers_prefix 'ConnectingGM_MoriGroups_DN.pdb'] ' &>' SAlogfile '&']);
        else

    end
pause(10);
end

end

return
%%%%%%%%%%% RUN AFTER ALL TRUESA and VALIDATION JOBS DONE -- this step
%%%%%%%%%%% is not required to continue on with full  brain BM
%%%%%%%%%%% normalization. This step is required for quality checks though.
%%%%%%%%%%%
%Tranform ind to mat
for s=1:size(subjectID, 2)
    cd([project_folder subjectID{s}]);
    dt6File=fullfile('dti06trilinrt', 'dt6.mat');

    %Get the pdb back to mat, incorporating the indices
    fg=dtiLoadFiberGroup(moriFile);
    try
    fid=fopen(fullfile(fileparts(dt6File), 'fibers', [allfibers_prefix  'ConnectingGM_MoriGroups_DN.ind']));
    fprintf( subjectID{s});
    %TrueSA solution are indices of fibers to keep--in the space of allConnectingGM_MoriGroups. PLUS1!!! (those indices count from 0)
    DN_ind=textscan(fid, '%d');DN_ind=DN_ind{1}; fclose(fid);
    dtiWriteFibersSubset(moriFile, DN_ind+1, fullfile(fileparts(dt6File), 'fibers', [allfibers_prefix  'ConnectingGM_MoriGroups_DN.mat']), 'allConnectingGM_MoriGroups_DN', fg.subgroup(DN_ind+1), fg.subgroupNames);

    catch
        fprintf('Not converged yet');
    end
end

return

%After the initial solution has been obtained, use its model.pid to run
%blue matter on the full set arranged as MoriTracts-then-otherTracts.
a=1;
while a
try
    a=(size(findstr('allConnectingGM_DN.ind', ls([project_folder filesep '*' filesep 'dti06trilinrt' filesep 'fibers' filesep allfibers_prefix  'ConnectingGM_DN.ind'])), 2)<112);
catch
    %in no *.ind yet
    a=1;
end
for s=1:size(subjectID, 2)
    %check that the initial guess job at least finished. system(grep) returns "0" if
    %there is a match found---oukey. (system error codes of 1 or 2 otherwise).. And also see if maybe this job is
    %already running...
        display([subjectID{s}]);

    [status, result]=system('ps -U elenary -u elenary |grep ''trueSA'' |wc|awk ''{print $1}''');
    if(str2num(result)<=(numProcessorsAvail+2)&&~system(['grep "converged" /azure/scr1/bluematterscratch/' subjectID{s} 'dti06trilinrt_' allfibers_prefix 'ConnectingGM_MoriGroups_DN' filesep '*true*']) && ~exist(['/azure/scr1/bluematterscratch/' subjectID{s} 'dti06trilinrt_' allfibers_prefix 'ConnectingGM_DN' filesep 'data_0.SBfloat'], 'file') && (system(['grep "converged" /azure/scr1/bluematterscratch/' subjectID{s} 'dti06trilinrt_' allfibers_prefix 'ConnectingGM_DN' filesep '*true*'])~=0))

        display(['Submitting all fibers Blue Matter for ' subjectID{s}]);
        cd([project_folder subjectID{s}]);
        initialGuess=['/azure/scr1/bluematterscratch/' subjectID{s} 'dti06trilinrt_' allfibers_prefix 'ConnectingGM_MoriGroups_DN' filesep 'pids' filesep 'model.pid'];
        system(['/biac3/wandell4/users/elenary/density_normalization/scripts/runSAmpi.sh ' project_folder filesep subjectID{s} filesep 'dti06trilinrt  ' [allfibers_prefix 'ConnectingGM.pdb'] ' ' [allfibers_prefix 'ConnectingGM_DN.pdb '] initialGuess ' &']);
    end
pause(50)
end

end

return

%Final solution -- Tranform ind to mat.  allConnectingGM_DN.mat: Mori Labels for the indices correspoinding to
%first N fibers can be obtained from allConnectingGM_MoriGroups, the rest
%should be "0";
for s=1:size(subjectID, 2)
    cd([project_folder subjectID{s}]);
    dt6File=fullfile('dti06trilinrt', 'dt6.mat');

    %Get the pdb back to mat, incorporating the indices
    fgWholeBrain=dtiLoadFibergroup(wholebrainGM);
    try
    fid=fopen(fullfile(fileparts(dt6File), 'fibers', [allfibers_prefix  'ConnectingGM_DN.ind']));
    fprintf( subjectID{s});
    %TrueSA solution are indices of fibers to keep--in the space of allConnectingGM. PLUS1!!! (those indices count from 0)
    DN_ind=textscan(fid, '%d');DN_ind=DN_ind{1}; fclose(fid);
    dtiWriteFibersSubset(wholebrainGM, DN_ind+1, fullfile(fileparts(dt6File), 'fibers', [allfibers_prefix  'ConnectingGM_DN.mat']), 'allConnectingGM_DN', fgWholeBrain.subgroup(DN_ind+1), fgWholeBrain.subgroupNames);

    catch
        fprintf('Not converged yet');
    end
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function ft_compile_mex(force)

% FT_COMPILE_MEX can be used for compiling most of the FieldTrip MEX files Note that
% this function does not put the MEX files in the correct location in the private
% folders, this is managed by a Bash script. In case you are not working with Git and
% you want to recompile the mex files for your platform, you can find all mex files
% for your platform and move them to a backup directory that is not on your MATLAB
% path. Subsequently you can rtun this function to recompile it on your platform with
% your compiler settings
%
% The standards procedure for compiling mex files is detailled on
% http://www.fieldtriptoolbox.org/development/guidelines/code#compiling_mex_files
%
% Please note that this script does NOT set up your MEX environment for you, so in
% case you haven't selected the C compiler on Windows yet, you need to type 'mex
% -setup' first to choose either the LCC, Borland or Microsoft compiler. If you want
% to use MinGW, you also need to install Gnumex (http://gnumex.sourceforget.net),
% which comes with its own procedure for setting up the MEX environment.
%
% The logic in this script is to first build a list of files that actually need
% compilation for the particular platform that MATLAB is running on, and then to go
% through that list. Functions are added to the list by giving their destination
% directory and (relative to that) the name of the source file (without the .c).
% Optionally, you can specify a list of platform this file needs to be compiled on
% only, and a list of platforms where you don't compile it on. Finally, you can give
% extra arguments to the MEX command, e.g., for including other c-sources or giving
% compiler flags.
%
% See also MEX

% Copyright (C) 2010, Stefan Klanke
%
% $Id$

if nargin<1
  force=false;
end

% Possible COMPUTER types
% GLNX86
% GLNXA64
% PCWIN
% PCWIN64
% MAC
% MACI
% MACI64

[ftver, ftpath] = ft_version;

L = [];
L = add_mex_source(L,'compat/matlablt2010b/@uint64','abs');
L = add_mex_source(L,'compat/matlablt2010b/@uint64','min');
L = add_mex_source(L,'compat/matlablt2010b/@uint64','max');
L = add_mex_source(L,'compat/matlablt2010b/@uint64','plus');
L = add_mex_source(L,'compat/matlablt2010b/@uint64','minus');
L = add_mex_source(L,'compat/matlablt2010b/@uint64','times');
L = add_mex_source(L,'compat/matlablt2010b/@uint64','rdivide');

L = add_mex_source(L,'@config/private','deepcopy');
L = add_mex_source(L,'@config/private','increment');
L = add_mex_source(L,'@config/private','setzero');

L = add_mex_source(L,'realtime/online_mri','ft_omri_smooth_volume');
L = add_mex_source(L,'realtime/src/acquisition/siemens/src', 'sap2matlab', [], [], 'siemensap.c -I../include');

L = add_mex_source(L,'src','ft_getopt');
L = add_mex_source(L,'src','read_16bit');
L = add_mex_source(L,'src','read_24bit');
L = add_mex_source(L,'src','read_ctf_shm', {'GLNX86'});  % only compile on GLNX86
L = add_mex_source(L,'src','write_ctf_shm', {'GLNX86'}); % only compile on GLNX86

L = add_mex_source(L,'src','lmoutr',[],[],'geometry.c -I.');
L = add_mex_source(L,'src','ltrisect',[],[],'geometry.c -I.');
L = add_mex_source(L,'src','plinproj',[],[],'geometry.c -I.');
L = add_mex_source(L,'src','ptriproj',[],[],'geometry.c -I.');
L = add_mex_source(L,'src','routlm',[],[],'geometry.c -I.');
L = add_mex_source(L,'src','solid_angle',[],[],'geometry.c -I.');
L = add_mex_source(L,'src','rfbevent',[],{'PCWIN', 'PCWIN64'},'d3des.c -I.'); % do not compile on WIN32 and WIN64
L = add_mex_source(L,'src','meg_leadfield1');
L = add_mex_source(L,'src','splint_gh');
L = add_mex_source(L,'src','plgndr');

L = add_mex_source(L,'src','ft_spike_sub_crossx');
L = add_mex_source(L,'src','rename');
L = add_mex_source(L,'src','getpid');

L = add_mex_source(L,'src','nanmean');
L = add_mex_source(L,'src','nanstd');
L = add_mex_source(L,'src','nanvar');
L = add_mex_source(L,'src','nansum');
L = add_mex_source(L,'src','nanstd');
L = add_mex_source(L,'src','det2x2');
L = add_mex_source(L,'src','inv2x2');
L = add_mex_source(L,'src','mtimes2x2');
L = add_mex_source(L,'src','sandwich2x2');
L = add_mex_source(L,'src','combineClusters');

% this one is located elsewhere
L = add_mex_source(L,'external/fileexchange','CalcMD5',[],[],'CFLAGS=''-std=c99 -fPIC''');

% this one depends on the MATLAB version
if ft_platform_supports('libmx_c_interface')
  % use the C interface
  L = add_mex_source(L,'src','mxSerialize_c');
  L = add_mex_source(L,'src','mxDeserialize_c');
else
  % use the C++ interface
  L = add_mex_source(L,'src','mxSerialize_cpp');
  L = add_mex_source(L,'src','mxDeserialize_cpp');
end

oldDir = pwd;
try
  compile_mex_list(L, ftpath, force);
catch
  % the "catch me" syntax is broken on MATLAB74, this fixes it
  me = lasterror;
  cd(oldDir);
  rethrow(me);
end
cd(oldDir);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%
% Use as
%   list = add_mex_source(list, directory, relName, includePlatform, excludePlatform, extras)
%
% The list is a structure array of directory names, source file names, and
% extra arguments required for the compilation of MEX files. This function will
% create a new element of this structure and append it to L.
%
%   directory       = target directory of the mex-file
%   relName         = source file relative to 'directory'
%   includePlatform   = list of platforms this MEX file should only be compiled for.
%                     use an empty matrix [] to compile for all platforms
%   excludePlatform = list of platforms this MEX file should NOT be compiled for.
%   extras          = extra arguments to the MEX command, e.g. additional source files
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function L = add_mex_source(L, directory, relName, includePlatform, excludePlatform, extras)

% Check if this file only needs compilation on certain platforms (including this one)
if nargin>3 && ~isempty(includePlatform)
  ok = false;
  for k=1:numel(includePlatform)
    if strcmp(includePlatform{k}, computer)
      ok = true;
      break;
    end
  end
  if ~ok
    return
  end
end

% Check if this file cannot be compiled on certain platforms (including this one)
if nargin>4 && ~isempty(excludePlatform)
  ok = true;
  for k=1:numel(excludePlatform)
    if strcmp(excludePlatform{k}, computer)
      return;
    end
  end
end

L(end+1).dir   = directory;
L(end).relName = relName;
if nargin>5
  L(end).extras = extras;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%
% Use as
%  compile_mex_list(L, baseDir)
%
% Compile a list of MEX files as determined by the input argument L.
% The second argument 'baseDir' is the common base directory for the
% files listed in L. The third argument is a flag that determines
% whether to force (re-)compilation even if the MEX file is up-to-date.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function compile_mex_list(L, baseDir, force)

for i=1:length(L)
  [relDir, name] = fileparts(L(i).relName);
  sfname1 = fullfile(baseDir, L(i).dir, [L(i).relName '.c']);
  sfname2 = fullfile(baseDir, L(i).dir, [L(i).relName '.cpp']);
  if exist(sfname1, 'file')
    sfname = sfname1;
    L(i).ext = 'c';
  elseif exist(sfname2, 'file')
    sfname = sfname2;
    L(i).ext = 'cpp';
  else
    sfname = '';
    L(i).ext = '';
    fprintf(1,'Error: source file for %s cannot be found.\n', L(i).relName);
    continue;
  end
  SF = dir(sfname);

  if ~force
    mfname = [baseDir filesep L(i).dir filesep name '.' mexext];
    MF = dir(mfname);
    if numel(MF)==1 && datenum(SF.date) <= datenum(MF.date)
      fprintf(1,'Skipping up-to-date MEX file %s/%s\n', L(i).dir, name);
      continue;
    end
  end
  fprintf(1,'Compiling MEX file %s/%s ...\n', L(i).dir, name);
  cd([baseDir '/' L(i).dir]);
  cmd = sprintf('mex %s.%s %s', L(i).relName, L(i).ext, L(i).extras);
  eval(cmd);
end


text:
function [ t_stop, y_stop ] = p06_stop ( neqn )

%*****************************************************************************80
%
%% P06_STOP returns the stopping point for problem p06.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real T_STOP, Y_STOP(NEQN), the final data.
%
  y_stop = zeros ( neqn, 1 );

  t_stop = 20.0;
  y_stop = [ ...
    0.6761876008576667; ...
    0.1860816099640036 ];

  return
end


text:
function qamsys(bin,f)

disp('========================================');
disp(' HAM DIEU CHE DICH BIEN & PHA: QAM');
disp(' VI DU: qamsys([0 0 0 1 1 1 0 1 1 1 1 1],3)');
disp('                Writer: Minh Dai Ca');
disp('========================================');

bin=[0 0 0 1 1 1 0 1 1 1 1 1];f=5;
t=0:2*pi/149:2*pi;

L = length(bin);bit1=ones(1,50);bit0=zeros(1,50);mbit=[];mcw=[];
sig000=cos(f*t);sig001=cos(f*t+pi/2);sig010=cos(f*t+pi);sig011=cos(f*t+3*pi/2);
sig100=2*sig000;sig101=2*sig001;sig110=2*sig010;sig111=2*sig011;

if 3*fix(L/3)~=L
    error('DO DAI CUA CHUOI bin PHAI LA BOI SO CUA 3');
end

for n=1:3:L;
    if bin(n)==0 && bin(n+1)==0 && bin(n+2)==0;
       cw=sig000;bit=[bit0 bit0 bit0];
    elseif bin(n)==0 && bin(n+1)==0 && bin(n+2)==1;
       cw=sig001;bit=[bit0 bit0 bit1];
    elseif bin(n)==0 && bin(n+1)==1 && bin(n+2)==0;
       cw=sig010;bit=[bit0 bit0 bit0];
    elseif bin(n)==0 && bin(n+1)==1 && bin(n+2)==1;
       cw=sig011;bit=[bit0 bit1 bit1];

    elseif bin(n)==1 && bin(n+1)==0 && bin(n+2)==0;
       cw=sig100;bit=[bit0 bit0 bit0];
    elseif bin(n)==1 && bin(n+1)==0 && bin(n+2)==1;
       cw=sig101;bit=[bit0 bit0 bit1];
    elseif bin(n)==1 && bin(n+1)==1 && bin(n+2)==0;
       cw=sig110;bit=[bit0 bit0 bit0];
    elseif bin(n)==1 && bin(n+1)==1 && bin(n+2)==1;
       cw=sig111;bit=[bit0 bit1 bit1];
    end
    mbit=[mbit bit];
   mcw=[mcw cw];
end
qam=mcw;
deqam=mcw;

subplot(3,1,1);plot(mbit,'r','linewidth',2);axis([0  50*L -0.5 1.5]);grid on;title('Data in');
subplot(3,1,2);plot(qam,'m','linewidth',1.5);axis([0  50*L -2.5 2.5]);grid on;title('PSK modulation');
subplot(3,1,3);plot(deqam,'g','linewidth',1.5);axis([0  50*L -2.5 2.5]);grid on;title('PSK demodulation,Data out');


text:
function dumpIfileHeader(ifileName)
% dumpIfileHeader(ifileName)
%
% Displays some of the GE header fields for the given I-file
% (eg. 'I.001') in human-readable format.
%
% REQUIRES:
%  GE2SPM tools- eg. addpath /usr/local/matlab/toolbox/mri/filters/GE2SPM/
%   -or-
%  matlab's DICOM tools (part of the image processing toolbox in matlab
%  version >= 6.5)
%
% HISTORY:
%  2002.06.07 RFD (bob@white.stanford.edu) wrote it.
%  2004.07.01 RFD: we now use Junjie's readImageHeader, so that this
%  function will work with either GE I-files or DICOM files.

if(nargin<1)
    help(mfilename);
    return;
end

[su_hdr,ex_hdr,se_hdr,im_hdr] = readImageHeader(ifileName);

disp(sprintf('Exam #: %d', ex_hdr.ex_no));
% GE timestamp is a unix timestamp- # of seconds since 1970.01.01 00:00:00
% eg. yr = round(1970+ex_hdr.ex_datetime/60/60/24/365.25); ...
disp(sprintf('Exam timestamp: %d', ex_hdr.ex_datetime));
disp(sprintf('Hospital name: %s', char(ex_hdr.hospname')));
disp(sprintf('Field strength: %0.1f Tesla', ex_hdr.magstrength/10000));
disp(sprintf('Patient ID: %s', char(ex_hdr.patid')));
disp(sprintf('Patient name: %s', char(ex_hdr.patname')));
disp(sprintf('Patient weight: %0.1f Kg', ex_hdr.patweight/1000));
disp(sprintf('Patient age: %0.0f years', ex_hdr.patage));
disp(sprintf('History: %s', char(ex_hdr.hist')));
disp(sprintf('Exam description: %s', char(ex_hdr.ex_desc')));

disp(sprintf('PSD iname: %s', char(im_hdr.psd_iname')));
disp(sprintf('Scan Time: %0.1f minutes', im_hdr.sctime/1000000/60));
disp(sprintf('FOV: %0.1f', im_hdr.dfov));
disp(sprintf('X,Y dim: %d, %d', im_hdr.dim_X, im_hdr.dim_Y));
disp(sprintf('X,Y pixels: %d, %d', im_hdr.imatrix_X, im_hdr.imatrix_Y));
disp(sprintf('nSlices: %d @ %0.4f mm thickness', im_hdr.slquant, im_hdr.slthick));
disp(sprintf('in-plane voxel size: x=%0.4f mm, y=%0.4fmm', im_hdr.pixsize_X, im_hdr.pixsize_Y));
disp(sprintf('TR: %0.1f ms', im_hdr.tr));
disp(sprintf('TI: %0.1f ms', im_hdr.ti));
disp(sprintf('TE: %0.1f ms (%d echos)', im_hdr.te, im_hdr.numecho));
disp(sprintf('flip: %d deg', im_hdr.mr_flip));
disp(sprintf('NEX: %f', im_hdr.nex));
disp(sprintf('SAR average: %f, peak: %f', im_hdr.saravg, im_hdr.sarpeak));

return;


text:
function cvx_optpnt = nonnegative( sx ) %#ok

%NONNEGATIVE   The nonnegative orthant.
%   NONNEGATIVE(SX), where SX is a valid size vector, creates an array
%   of size SX and constrains each element to be nonnegative. Therefore,
%   given the declaration
%      variable x(sx)
%   the constraint
%      x == nonnegative(sx);
%   is equivalent to
%      x >= 0;
%   Obviously, the inequality form is simpler and is preferred in most
%   circumstances.
%
%   Disciplined convex programming information:
%       NONNEGATIVE is a cvx set specification. See the user guide for
%       details on how to use sets.

narginchk(1,1);

[ temp, sx ] = cvx_check_dimlist( sx, true );
if ~temp,
    error( 'Argument must be a non-empty dimension vector.' );
end

cvx_begin set
    variables x( sx )
    if all( sx ~= 0 ),
        [ tx, dummy ] = find( cvx_basis( x ) ); %#ok
        newnonl( cvx_problem, 'nonnegative', tx(:) );
    end
cvx_end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
% sgwt_setpath : Set paths for SGWT toolbox
%
% user should set SGWT_ROOT to be directory where sgwt_toolbox is installed

% This file is part of the SGWT toolbox (Spectral Graph Wavelet Transform toolbox)
% Copyright (C) 2010, David K. Hammond.
%
% The SGWT toolbox is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% The SGWT toolbox is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with the SGWT toolbox.  If not, see <http://www.gnu.org/licenses/>.

SGWT_ROOT=[fileparts(mfilename('fullpath')),filesep];

sgwt_relpathlist={'',...
                  'utils',...
                  'demo',...
                  'mex',...
                 };
fprintf('Welcome to sgwt_toolbox. SGWT root directory is %s\n',SGWT_ROOT);
for k=1:numel(sgwt_relpathlist)
  sgwt_tmp_pathname=[SGWT_ROOT,sgwt_relpathlist{k}];
  fprintf('adding path %s\n',sgwt_tmp_pathname);
  addpath(sgwt_tmp_pathname);

end
clear sgwt_relpathlist sgwt_tmp_pathname

repo: gspbox


text:
function [V,F] = bwstencil(im,varargin)
  % BWSTENCIL Create a 3D-printable stencil from an input binary image.
  %
  % [V,F] = bwstencil(im)
  % [V,F] = bwstencil(im,'ParameterName',ParameterValue, ...)
  %
  % Inputs:
  %   im  h by w by 1 binary image
  %   Optional inputs (see bwmesh.m)
  % Outputs:
  %   V  #V by 3 list of vertex positions with V(:,3) = 0 or 1 for top or
  %   bottom layer
  %
  % Example:
  %   % Read in image, convert to binary, shrink, invert
  %   im = ~imresize(im2bw(imread('hans-hass.jpg')),0.25);
  %   % Pad image to create boundary frame
  %   im = padarray(im,ceil(max(size(im))*0.125*[1;1]));
  %   % Generate stencil of thickness = 1px
  %   [V,F] = bwstencil(im,'Tol',1.5,'SmoothingIters',50);
  %   % Change thickness to 2.5% size of image
  %   V(:,3) = V(:,3) * 0.025*max(size(im));
  %   % render result:
  %   subplot(2,1,1);
  %   imshow(im);
  %   subplot(2,1,2);
  %   tsurf(F,V,'EdgeColor','none','FaceColor',[0.3 0.4 1.0],'FaceAlpha',0.9);
  %   l = light('Style','infinite','Position',[0.3 0.2 1]);
  %   axis equal;
  %
  % See also: bwmesh
  %
  if ~islogical(im)
    warning('Converting non-binary input to binary');
    im = im2bw(im);
  end
  imh = imfill(im,'holes');
  if any(im(:) ~= imh(:))
    warning('input contained islands, removing...');
    im = imh;
  end
  im = imresize(im,2,'nearest');
  [V,F] = bwmesh(~im,varargin{:});
  [V,F] = extrude(V,F);
  [V,F] = remesh_planar_patches(V,F);
end

author: alecjacobson

repo: gptoolbox


text:
function O = GetOPTI(B,opts)
%Build OPTI Object from SymBuilder Object
%
%   Called By SYMBUILDER Class

%   Copyright (C) 2012-2013 Jonathan Currie (IPL)

if(nargin < 2 || isempty(opts)), opts = symbset; else, opts = symbset(opts); end

%Check object is built
if(~IsBuilt(B))
    error('Please build the object using Draft() or Build()');
end

%Check we have an objective (no good without!)
if(B.noObjs == 0)
    error('You can only construct an OPTI object when you have specified an objective!');
elseif(B.noObjs > 1)
    error('You can only construct an OPTI object with a single objective');
end

%Clear existing callback function from memory
clear symb_cb symb_ccb

%Separate objective and constraint indices
olin = B.objlin(B.indobj);
clin = B.conlin; clin(B.indobj) = [];

%Default options
if(~isempty(opts.optiOpts))
    Oopts = optiset(opts.optiOpts,'warnings','off'); %overrides solver and display settings
    %check if user has specified a solver/display
    if(~strcmpi(opts.solver,'auto'))
        optiwarn('OPTI:SymBSolver','As you have supplied optiOpts to symbset, your solver selection within symbset has been overwritten by the solver specified in optiOpts.');
    end
    if(~strcmpi(opts.display,'off'))
        optiwarn('OPTI:SymBDisplay','As you have supplied optiOpts to symbset, your display setting within symbset has been overwritten by the display setting specified in optiOpts.');
    end
else
    Oopts = optiset(opts.optiOpts,'warnings','off','solver',opts.solver,'display',opts.display);
end

%Check for Linear Problem
if(olin == 1 && all(clin == 1))
    %Build OPTI object
    O = opti(GetLinProb(B),'options',Oopts);
%Check for Quadratic Pprogram
elseif(olin == 2 && all(clin == 1))
    %Build OPTI object
    O = opti(GetQuadProb(B),'options',Oopts);
%Check for QCQP/QCLP
elseif(olin <= 2 && all(clin <= 2))
    %Build OPTI object
    O = opti(GetQuadProb(B),'options',Oopts);
%Check for NLP
elseif(olin <= 3 && all(clin <= 3))
    %If we are building an MINLP, add BONMIN linear information
    if(any(B.xtype ~= 'C'))
        %Constraint Linearity
        cons_lin = clin;
        cons_lin(clin > 1) = 0;
        %Variable linearity
        varlin = ones(size(B.jac,2),1);
        nlvars = symvar(B.jac); %Variables in the jacobian are nonlinear WHAT ABOUT GRADIENT?
        for i = 1:length(nlvars)
            varlin(logical(nlvars(i) == B.vars)) = 0;
        end
        if(~isempty(Oopts.solverOpts))
            Oopts = optiset(Oopts,'warnings','off','solverOpts',bonminset(Oopts.solverOpts,'cons_lin',cons_lin,'var_lin',varlin));
        else
            Oopts = optiset(Oopts,'warnings','off','solverOpts',bonminset('cons_lin',cons_lin,'var_lin',varlin));
        end
    end
    %Build OPTI object
    O = opti(GetNLProb(B,opts),'options',Oopts);
else
    error('problem type not yet supported');
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function cprops = edl_check_internalindices(props)
%EDL_CHECK_INTERNALINDICES Checks the consistency of internal indices
%
% $ Syntax $
%   - cprops = edl_check_internalindices(props)
%
% $ Arguments $
%   - props:        the struct array of property entries
%   - cprops:       the converted array with all indices coverted to
%                   numeric
%
% $ History $
%   - Created by Dahua Lin, on Aug 14, 2006
%

cprops = props;

if ~isempty(props)

    n = length(props);

    if ~isfield(props, 'internal_index')
        error('edl:parseerror', ...
            'The entries do not have the required field: internal_index');
    end

    for i = 1 : n
        curidx = str2double(props(i).internal_index);
        if curidx ~= i
            error('edl:parseerror', ...
                'Internal index inconsistency');
        end
        cprops(i).internal_index = curidx;
    end

end

author: lmthang


text:
function [t, x, u] = fst_mpc(fst, output_data, varargin)
    warning off all;
%%
    fst.net_load = fst.load *2 - fst.PV/3 - fst.wind; %Power Balance

    %   Solve the optimal control problem
    t_Start = tic;
    [fst.u0, V_current, exitflag, output] = solveOptimalControlProblem( fst );
    t_Elapsed = toc( t_Start );

    if ( fst.iprint >= 1 )
        printSolution(fst.printClosedloopData, fst.mpciter, fst.xmeasure, fst.u0, fst.iprint, exitflag, t_Elapsed,output_data);
    end

    % Step (3):
    %Save current information:
    t = output;
    x = computeOpenloopSolution(fst,fst.u0);
    u = fst.u0;

end

author: juchengquan


text:
% EASY12   Detailed explanations of a Matlab implementation
%          of the Lambda method. The notation follows the one
%          introduced in Strang and Borre, pages 495--499

%Kai Borre 30-06-2008
%Copyright (c) by Kai Borre
%$Revision: 1.0 $  $Date: 2008/06/30  $

fprintf('\n')
echo on
% Test example
hat_I = [5.45;3.1;2.97];
%hat_I =[0; 1];
n = size(hat_I,1);
Q_hat_I = [6.29 5.978 .544; 5.978 6.292 2.34; .544 2.34 6.288];
%Q_hat_I = [53.40 38.40;38.40 28.00];
echo off

% Given float ambiguities
fprintf('\nFloat ambiguities hat_I')
for i = 1:n
    fprintf('\n%12.3f', hat_I(i,1))
end

% Original covariance matrix for hat_I
fprintf('\n\nCovariance matrix for hat_I')
for i = 1:n
    fprintf('\n')
    for j = 1:n
        fprintf('%12.3f', Q_hat_I(i,j))
    end
end

fprintf('\n\n')
echo on
% Shift of hat_I in order to secure that -1 < hat_I <= +1
echo off

shifts = hat_I-rem(hat_I,1);
fprintf('\n\nInteger shifts of hat_I')
for i = 1:n
    fprintf('\n%12.0f', shifts(i,1))
end

% Remainders of hat_I
hat_I_r = rem(hat_I,1);
fprintf('\n\nRemainders of hat_I')
for i = 1:n
    fprintf('\n%12.3f', hat_I_r(i,1))
end

% L^T*D*L factorization of Q_hat_I
% For consistency we do not use the Matlab ldl function
[L,D] = ldldecom(Q_hat_I);
fprintf('\n\nQ_hat_I is factorized into L^T*D*L')
fprintf('\n\nThe lower triangular L')
for i = 1:n
    fprintf('\n')
    for j = 1:n
        fprintf('%12.3f',L(i,j))
    end
end
fprintf('\n\nThe diagonal matrix D\n')
for i = 1:n
    fprintf('%12.3f',D(i))
end

% Computing the size of the search volume
chi2 = chistart(D,L,hat_I_r,1);
fprintf('\n\nchi^2 = %5.3f\n',chi2)

% Doing the decorrelation
[Q_bar_I,Z,L_t,D_t] = decorrel(Q_hat_I,hat_I_r);
fprintf('\nInteger transformation matrix Z')
for i = 1:n
    fprintf('\n')
    for j = 1:n
        fprintf('%12.0f',Z(i,j))
    end
end

% hat_I transformed: I_s = Z'*hat_I_r
fprintf('\n\nTransformed, shifted ambiguities I_s = Z^T*hat_I_r')
I_s = Z'*hat_I_r;
for i = 1:n
    fprintf('\n%12.3f',I_s(i))
end

% The transformed, decorrelated covariance matrix: Q_bar_I = Z^T*Q_hat_I*Z
fprintf('\n\nThe transformed, decorrelated covariance matrix Q_bar_I = Z^T*Q_hat_I*Z.')
fprintf('\n\nQ_bar_I = Z^T*Q_hat_I*Z')
for i = 1:n
    fprintf('\n')
    for j = 1:n
        fprintf('%12.3f',Q_bar_I(i,j))
    end
end
fprintf('\nNote the deminished off-diagonal terms!')
fprintf('\n\nThe lower triangular L_t used in the search')
for i = 1:n
    fprintf('\n')
    for j = 1:n
        fprintf('%12.3f',L_t(i,j))
    end
end
fprintf('\n\nThe diagonal matrix D_t used in the seach\n')
for i = 1:n
    fprintf('%12.3f',D_t(i))
end

fprintf('\n\n')
echo on
% Determining the size of the search volume for the transformed L_t and D_t
echo off

chi2 = chistart(D_t,L_t,I_s,1);
fprintf('\n\nchi^2 for the transformed problem  =  %5.3f\n',chi2)

fprintf('\nThe search domain is defined as')
fprintf('\n(I-hat_I_r)^T*(Z^T*Q_hat_I*Z)*(I-hat_I_r) < chi^2,  for I integer')

[bar_I,sqnorm,ierr] = lsearch(I_s,L_t,D_t,chi2,1);
%bar_I = (bar_I' * inv(Z))'+shifts;
bar_I = (bar_I'/Z)' + shifts;
fprintf('\n\nFixed ambiguities bar_I')
for i = 1:n
    fprintf('\n%12.0f', bar_I(i,1))
end
fprintf('\n\n')
%%%%%%%%%%%%%%%%%%%%%%%%%%% end easy12.m  %%%%%%%%%%%%

author: yandld


text:
function model = knReg(X, t, lambda, kn)
% Gaussian process (kernel) regression
% Input:
%   X: d x n data
%   t: 1 x n response
%   lambda: regularization parameter
% Output:
%   model: trained model structure
% Written by Mo Chen (sth4nth@gmail.com).
if nargin < 4
    kn = @knGauss;
end
if nargin < 3
    lambda = 1e-2;
end
K = knCenter(kn,X);
tbar = mean(t);
U = chol(K+lambda*eye(size(X,2)));    % 6.62
a = U\(U'\(t(:)-tbar));               % 6.68

model.kn = kn;
model.a = a;
model.X = X;
model.tbar = tbar;
%% for probability prediction
y = a'*K+tbar;
beta = 1/mean((t-y).^2);              % 3.21
alpha = lambda*beta;           % lambda=a/b P.153 3.55
model.alpha = alpha;
model.beta = beta;
model.U = U;

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function [VCV,A,B,scores,hess,gross_scores]=robustvcv(fun,theta,nw,varargin)
% Compute Robust Variance Covariance matrix numerically, including
% Newey-West style score covariance using 2-sided derivatives
%
% USAGE:
%     [VCV,A,B,SCORES,HESS,GROSS_SCORES]=robustvcv(FUN,THETA,NW,VARARGIN)
%
% INPUTS:
%     FUN           - Function name ('fun') or function handle (@fun) which will
%                       return the sum of the log-likelihood (scalar) as the 1st output and the individual
%                       log likelihoods (T by 1 vector) as the second output.
%     THETA         - Parameter estimates at the optimum, usually from fmin*
%     NW            - Number of lags to consider in Newey-West covariance.
%                       Normally set to 0
%     VARARGIN      - Other inputs to the log-likelihood function, such as data
%
% OUTPUTS:
%     VCV           - Estimated robust covariance matrix (see White 1994)
%     A             - A portion of robust covariance
%     B             - B portion of robust covariance
%     SCORES        - T x num_parameters matrix of scores
%     HESS          - Estimated Hessian (Expectation of second derivative)
%     GROSS_SCORES  - Numerical scores (1 by num_parameters) of the objective function, usually for diagnostics
%
% COMMENTS:
%     This function simplifies calculating sandwich covariance estimators for (Q)MLE estimation

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 1    Date: 9/1/2005

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Argument Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if size(theta,1)<size(theta,2)
    theta=theta';
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Argument Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

k=length(theta);
h=max(abs(theta*eps^(1/3)),1e-8);
h=diag(h);

[~,like]=feval(fun,theta,varargin{:});

t=length(like);

LLFp=zeros(k,1);
LLFm=zeros(k,1);
likep=zeros(t,k);
likem=zeros(t,k);
for i=1:k
    thetaph=theta+h(:,i);
    [LLFp(i),likep(:,i)]=feval(fun,thetaph,varargin{:});
    thetamh=theta-h(:,i);
    [LLFm(i),likem(:,i)]=feval(fun,thetamh,varargin{:});
end

scores=zeros(t,k);
gross_scores=zeros(k,1);
h=diag(h);
for i=1:k
    scores(:,i)=(likep(:,i)-likem(:,i))./(2*h(i));
    gross_scores(i)=(LLFp(i)-LLFm(i))./(2*h(i));
end

hess=hessian_2sided(fun,theta,varargin{:});
A=hess/t;
hess=A;
Ainv=A^(-1);
if nw==0
    % VCV=A^(-1)*B*A^(-1)/t;
    B=cov(scores);
    VCV=(Ainv*B*Ainv)/t;
else
    B=covnw(scores,nw);
    VCV=(Ainv*B*Ainv)/t;
end

author: bashtage

repo: mfe-toolbox


text:
%
% Removes graph nodes that are within one pixel of eachother
%
% Input
%   T: graph struct that is the output of the tracer
%
% Output
%   T:
%
function T = graph_merge(T)

    % make all merge moves
    tormv = true;
    while tormv
       [T,tormv] = remove_conncomp(T);
    end

end

% make one removal
function [M,tormv] = remove_conncomp(M)

    G = M.G;
    D = squareform(pdist(G));
    EE = D < 1.1; % make an edge between two nodes that are very close
    [S,C] = graphconncomp(sparse(EE),'Directed',false);
    % C indicates which component each node belongs to

    tormv= false;
    for i=1:S
       sel = C==i;
       ns = sum(sel);
       if ns > 1
          tormv = true;
          break
       end
    end

    % break when necessary
    if ~tormv, return; end

    %E_sel_to_sel = M.E(sel,sel);
    E_sel_to_not = M.E(sel,~sel);
    EI_sel_to_sel = M.EI(sel,sel);
    EI_sel_to_not = M.EI(sel,~sel);

    % create the new row in adjacency matrix, which
    % includes all the nodes connected to the selected set/pair
    new_row_E = any(E_sel_to_not,1);
    new_row_EI = cell(size(new_row_E));
    for i=1:size(EI_sel_to_not,2);
        new_row_EI{i} = tovec(EI_sel_to_not(:,i));
    end

    % list_ei_sel_to_not = tovec(EI_sel_to_not);
    list_ei_sel_to_sel = tovec(EI_sel_to_sel);
    list_ei_sel_to_sel = unique(list_ei_sel_to_sel);

    % see which of the edges, from sel to sel, should be removed.
    % we don't want to remove edges that make a long arc and then return
    nss = length(list_ei_sel_to_sel);
    rmv_edge = false(nss,1);
    for i=1:nss
       indx = list_ei_sel_to_sel(i);
       edge_length = size(M.S{indx},1);
       if edge_length == 2
          rmv_edge(i) = true;
       end
    end
    keep_ei_sel_to_sel =  list_ei_sel_to_sel(~rmv_edge);
    rmv_ei_sel_to_sel = list_ei_sel_to_sel(rmv_edge);

    % remove node positions
    G_sel = M.G(sel,:);
    new_node = mean(G_sel,1);
    Grmv = M.G(sel,:);
    M.G(sel,:) = [];
    M.G = [M.G; new_node];
    M.n = M.n - sum(sel) + 1;

    % modify adjacency matrix with the new row
    M.E(sel,:) = [];
    M.E(:,sel) = [];
    M.E(end+1,:) = new_row_E;
    M.E(:,end+1) = [new_row_E,false];
    if ~isempty(keep_ei_sel_to_sel)
        M.E(end,end) = true;
    end

    % modify the list of edges with each cell
    M.EI(sel,:) = [];
    M.EI(:,sel) = [];
    M.EI(end+1,:) = new_row_EI;
    M.EI(:,end+1) = [new_row_EI, {keep_ei_sel_to_sel} ];

    % modify the paths
    M.S(rmv_ei_sel_to_sel) = [];
    for i=1:numel(M.EI)
       mycell = M.EI{i};

       % update counts in adjacency matrix
       for j=1:length(mycell)
          el = mycell(j);
          M.EI{i}(j) = el - sum(rmv_ei_sel_to_sel < el);
       end
    end

    % update the paths, stored in M.S, such that all
    % nodes we have now replaced are updated with their
    % new coordinates
    for i=1:length(M.S)
        for j=1:size(Grmv,1);
            node_rmv = Grmv(j,:);
            d = pdist2(node_rmv,M.S{i});
            swap = d < .001;
            nswap = sum(swap);
            M.S{i}(swap,:) = repmat(new_node,[nswap 1]);
        end
    end

    assert_valid_graph(M);
end

author: brendenlake

repo: BPL

save_path: github-repos/MATLAB/brendenlake-BPL


text:
function [tfr,t,f] = tfrppage(x,t,N,h,trace);
%TFRPPAGE Pseudo Page time-frequency distribution.
%	[TFR,T,F]=TFRPPAGE(X,T,N,H,TRACE) computes the Pseudo Page
%	distribution of a discrete-time signal X, or the
%	cross Pseudo Page representation between two signals.
%
%	X     : signal if auto-PPage, or [X1,X2] if cross-PPage.
%	T     : time instant(s)          (default : 1:length(X)).
%	N     : number of frequency bins (default : length(X)).
%	H     : frequency smoothing window, H(0) being forced to 1.
%	                                 (default : Hamming(N/4)).
%	TRACE : if nonzero, the progression of the algorithm is shown
%	                                 (default : 0).
%	TFR   : time-frequency representation. When called without
%	        output arguments, TFRPPAGE runs TFRQVIEW.
%	F     : vector of normalized frequencies.
%
%	Example :
%	 sig=fmlin(128,0.1,0.4); tfrppage(sig);
%
%	See also all the time-frequency representations listed in
%	 the file CONTENTS (TFR*)

%	F. Auger, May-August 1994, July 1995.
%	Copyright (c) 1996 by CNRS (France).
%
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

[xrow,xcol] = size(x);
if (nargin < 1),
 error('At least 1 parameter is required');
elseif nargin<=2,
 N=xrow;
end;

hlength=floor(N/4);
if (rem(hlength,2)==0),
 hlength=hlength+1;
end;

if (nargin == 1),
 t=1:xrow; h = tftb_window(hlength); trace=0;
elseif (nargin == 2 | nargin == 3),
 h = tftb_window(hlength); trace = 0;
elseif (nargin == 4),
 trace = 0;
end;

if (N<0),
 error('N must be greater than zero');
end;
[trow,tcol] = size(t);
if (xcol==0)|(xcol>2),
 error('X must have one or two columns');
elseif (trow~=1),
 error('T must only have one row');
elseif (2^nextpow2(N)~=N),
 fprintf('For a faster computation, N should be a power of two\n');
end;

[hrow,hcol]=size(h); Lh=(hrow-1)/2; h=h/h(Lh+1);
if (hcol~=1)|(rem(hrow,2)==0),
 error('H must be a smoothing window with odd length');
end;

tfr= zeros (N,tcol) ;
if trace, disp('Pseudo Page distribution'); end;
for icol=1:tcol,
 ti= t(icol); tau=0:min([N-1,Lh,ti-1]);
 indices= rem(N+tau,N)+1;
 if trace, disprog(icol,tcol,10); end;
 tfr(indices,icol)=h(Lh+1+tau).*x(ti,1).*conj(x(ti-tau,xcol));
end;
if trace, fprintf('\n'); end;
tfr= real(fft(tfr));

if (nargout==0),
 tfrqview(tfr,x,t,'tfrppage',h);
elseif (nargout==3),
 if rem(N,2)==0,
  f=[0:N/2-1 -N/2:-1]'/N;
 else
  f=[0:(N-1)/2 -(N-1)/2:-1]'/N;
 end;
end;

author: HeLiangHIT


text:
startup

% specify your paths to the datasets
name = {'VV' 'LIME' 'NPE' 'NPE-ex1' 'NPE-ex2' 'NPE-ex3' 'MEF' 'DICM'};
dataset = strcat('data', filesep, name, filesep, '*.*');

% specify methods and metrics
method = {@multiscaleRetinex @dong @npe @lime @mf @srie @BIMEF};
metric = {@loe100x100};
% metric = {@loe100x100 @vif}; % NOTE matlabPyrTools is required to run VIF metric (vif.m).

for n = 1:numel(dataset); data = dataset{n};
    data,
    Test = TestImage(data);
    Test.Method = method;
    Test.Metric = metric;

    % run test and display results
    Test,

    % save test to a .csv file
    save(Test); % %save(Test, ['TestReport__' name{n} '.csv']);
end

author: baidut

repo: BIMEF

save_path: github-repos/MATLAB/baidut-BIMEF


text:
% -------------------------------------------------------------------------
% runduffsr.m
% Numerically solve the driven duffing oscillator with noise. Compute
% the averaged signal and noise amplitude spectra for varing noise
% strength. Plot signal-to-noise vs noise strength and time-series, phase
% space, well occupation distributions, and amplitude spectra for four
% chosen noise strengths.
% Dependencies: duff.m and noise.m
% -------------------------------------------------------------------------

clear all

tic % start timer

global p w

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

X0 = [0 0]'; % initial conditions
t0 = 0; % initial time
tf = 1000; % final time
dt = 0.05; % time step
T = (t0:dt:tf); % time vector
options = [];

f = (1:length(T))/(length(T)*dt); % frequency
f = f(1:(ceil(length(T)/2))); % first half of frequency axis

win = hann(length(T)); % hanning window

wmax = 10; % noise cut-off frequency
alpha = 0.5; % strength of non-linear term
beta = 1; % strength of linear term
gam = 0.5; % damping strength
A = 0.38; % driving amplitude
w0 = 0.09; % driving frequency

N = 1000; % number of terms in noise sum

d = 0.05:0.05:2; % noise strength

% noise strength values for plots
dc = [0.1 0.5 1 2];

dw0 = w0*0.25; % integration range (w0+/- 25%)

M = 4; % number of terms in average (>= 2)

% initialize arrays
x = zeros(length(T),length(dc));
v = zeros(length(T),length(dc));
lc = zeros(1,length(dc));

Fd = zeros((1+length(T))/2,M);
Fdn = zeros((1+length(T))/2,M);

Sd = zeros((1+length(T))/2,length(d));
Sdn = zeros((1+length(T))/2,length(d));

Pd = zeros(1,length(d));
Pdn = zeros(1,length(d));

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

for l = 1:length(d) % loop over noise strength array

n = d(l);

for j = 1:M % average amplitudes

p = unifrnd(0,2*pi,1,N); % random phases in noise
w = unifrnd(0,wmax,1,N); % random frequencies in noise

% numerically solve duffing oscillator
[t,y]=ode23(@duff,T,X0,options,gam,alpha,beta,A,w0,n);

% create x and v vectors at various noise strengths (d = 0.1,0.5,1.5,4)
% and get index for plots
if j == 1
    for i = 1:length(dc)
        if n == dc(i)
            x(:,i) = y(:,1);
            v(:,i) = y(:,2);
            lc(i) = l;
            break
        end
    end
end

% compute positive fft of windowed time-dependent signal and noise
% amplitudes
Fdtemp = fft(win.*y(:,1))/length(t); % signal fft
Fdntemp = fft(sqrt(n)*win.*noise(t)')/length(t); % noise fft

Fd(:,j) = Fdtemp(1:(ceil(length(t)/2))); % positive half of signal fft
Fdn(:,j) = Fdntemp(1:(ceil(length(t)/2))); % positive half of noise frequency

clear Fdtemp Fdntemp

end

% average frequency dependent signal and noise amplitude
Sd(:,l) = sum(abs(Fd'))/M; % averaged signal amplitude
Sdn(:,l) = sum(abs(Fdn'))/M; % averaged noise amplitude

i = 1;

% find signal and noise amplitude around driving frequency
for k = 1:length(f)
    if 2*pi*f(k) >= w0-dw0 && 2*pi*f(k) <= w0+dw0
        fr(i) = 2*pi*f(k); % frequency
        fd(i) = Sd(k,l); % signal amplitude
        fdn(i) = Sdn(k,l); % noise amplitude
        i = i + 1;
    end
end

% integrate signal and noise amplitude around driving frequency
Pd(l) = trapz(fr,fd); % signal "power"
Pdn(l) = trapz(fr,fdn); % noise "power"

clear fr fd fdn

end

% calculate and smooth snr (dB)
SNR = smooth(10*log(Pd./Pdn),0.1);

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

figure(1) % plot snr vs noise strength
plot(d,SNR,'-ob')
xlabel('d'); ylabel('SNR (dB)');

% plots for various noise strengths (d = 0.1,0.5,1.5,4)
figure(2) % plot time-series
subplot(4,1,1)
plot(t,x(:,1))
xlabel('t'); ylabel('x');
subplot(4,1,2)
plot(t,x(:,2))
xlabel('t'); ylabel('x');
subplot(4,1,3)
plot(t,x(:,3))
xlabel('t'); ylabel('x');
subplot(4,1,4)
plot(t,x(:,4))
xlabel('t'); ylabel('x');

figure(3) % plot phase-space
subplot(2,2,1)
plot(x(:,1),v(:,1))
xlabel('x'); ylabel('v');
subplot(2,2,2)
plot(x(:,2),v(:,2))
xlabel('x'); ylabel('v');
subplot(2,2,3)
plot(x(:,3),v(:,3))
xlabel('x'); ylabel('v');
subplot(2,2,4)
plot(x(:,4),v(:,4))
xlabel('x'); ylabel('v');

figure(4) % plot histogram of well distribution
subplot(2,2,1)
hist(x(:,1),100)
xlabel('x');
subplot(2,2,2)
hist(x(:,2),100)
xlabel('x');
subplot(2,2,3)
hist(x(:,3),100)
xlabel('x');
subplot(2,2,4)
hist(x(:,4),100)
xlabel('x');

figure(5) % plot noise subtracted amplitude spectra
plot(2*pi*f,Sd(:,lc(1))-Sdn(:,lc(1)),2*pi*f,Sd(:,lc(2))-Sdn(:,lc(2)),...
     2*pi*f,Sd(:,lc(3))-Sdn(:,lc(3)),2*pi*f,Sd(:,lc(4))-Sdn(:,lc(4)))
legend('d = 0.1','d = 0.6','d = 1.8','d = 2.4')
xlabel('\omega'); ylabel('|A|');

toc % end timer


text:
% ASTROTIK by Francesco Santilli
% R2BP (Restricted Two Bodies Problem)
% t2x computes orbital position vectors of a time span
%
% Usage: X = t2x(orbits,t0,t,mi)
%
% where: orbits(k,:) = [p e i o w] = 3d orbit elements
%        orbits(k,:) = [p e w]     = 2d orbit elements
%           p = semi-latus rectum [L] (p>0)
%           e = eccentricity [-] (e>=0)
%           i = inclinaion [rad]
%           o = raan [rad]
%           w = argument of perifocus [rad]
%        t0(k) = time of periapsis [T]
%        t(j) = time [T]
%        mi = gravitational parameter [L^3*T^-2] (mi>0)
%        X(j,:,k) = position [L]

function X = t2x(orbits,t0,t,mi)

    if ~(nargin == 4)
        error('Wrong number of input arguments.')
    end

    [K,D] = check(orbits,2);
    check(t0,1)
    J = check(t,1);
    check(mi,0)

    if D < 3
        error('Wrong size of input arguments.')
    end

    p = orbits(:,1);
    e = orbits(:,2);
    c = sqrt(mi./p.^3);
    f = t2f(t0,c,e,t);

    d3 = (D==5);
    X = zeros(J,2+d3,K);
    for k = 1:K
        X(:,:,k) = f2x(orbits(k,:),f(:,k));
    end

end


text:
function zList=genAllNestedZParenth(n,algorithm)
%%GENALLNESTEDZPARENTH0 Generate all nested sets of n parenthesis pairs,
%           recording only the indices of the left parantheses.
%
%INPUTS: n The number of nested paranthesis pairs; n>=0. Passing zero will
%          result in an empty matrix being returned.
% algorithm An optional parameter specifying the algorithm that should be
%          used to generate the parenthesis pairs. Possible values are:
%          0 (The default if omitted or an empty matrix is passed) Use
%            the algorithm of Problem 2 of Section 7.2.1.6 of [1].
%          1 Use algorithm N of Section 7.2.1.6 of [1].
%
%OUTPUTS: zList The nXnumPairs list of indices where the left parenthesis
%               goes.
%
%There is a total of CatalanNumber(n) possible arrangements of nested
%parentheses.
%
%EXAMPLE:
%To recreate the values in Table I of Section 7.2.1.6 of [1],
%one can run
% zList=genAllNestedZParenth(4)
%
%REFERENCES:
%[1] D. E. Knuth, The Art of Computer Programming. Vol. 4A: Combinatorial
%    Algorithms, Part I, Boston: Addison-Wesley, 2011.
%
%October 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2||isempty(algorithm))
    algorithm=0;
end

if(n==0)
    zList=[];
    return;
end

switch(algorithm)
    case 0
        zList=genAllNestedZParenth0(n);
    case 1
        zList=genAllNestedZParenth1(n);
    otherwise
        error('Unknown algorithm specified.')
end

end

function zList=genAllNestedZParenth0(n)
%%GENALLNESTEDZPARENTH0 Generate all nested sets of n parenthesis pairs,
%           recording only the indices of the left parantheses using the
%           algorithm in Problem 2 of Section 7.2.1.6 of [1].
%
%REFERENCES
%[1] D. E. Knuth, The Art of Computer Programming. Vol. 4A: Combinatorial
%    Algorithms, Part I, Boston: Addison-Wesley, 2011.
%
%October 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.

if(n==1)
    zList=1;
    return;
end

numBinTrees=CatalanNumber(n);
%Allocate space.
zList=zeros(n,numBinTrees);

%Step T1 Initialize
z=2*(0:n).'-1;

for curTree=1:numBinTrees
    %Step T2
    zList(:,curTree)=z(2:(n+1));

    %Step T3
    if(z(n-1+1)<z(n+1)-1)
        z(n+1)=z(n+1)-1;
        continue;
    end

    %Step T4, find j.
    j=n-1;
    z(n+1)=2*n-1;
    while(z(j-1+1)==z(j+1)-1)
        z(j+1)=2*j-1;
        j=j-1;
    end

    if(j==1)
        break;
    end
    z(j+1)=z(j+1)-1;
end

end

function zList=genAllNestedZParenth1(n)
%%GENALLNESTEDZPARENTH0 Generate all nested sets of n parenthesis pairs,
%           recording only the indices of the left parantheses using
%           Algorithm N of Section 7.2.1.6 of [1].
%
%REFERENCES
%[1] D. E. Knuth, The Art of Computer Programming. Vol. 4A: Combinatorial
%    Algorithms, Part I, Boston: Addison-Wesley, 2011.
%
%October 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.

numBinTrees=CatalanNumber(n);
%Allocate space.
zList=zeros(n,numBinTrees);

%Step N1
z=2*(1:n).'-1;
g=2*(1:n).'-2;

for curTree=1:numBinTrees
    %Step N2
    zList(:,curTree)=z;
    j=n;

    %Step N3, Find j.
    while(z(j)==g(j))
        g(j)=bitxor(g(j),1);
        j=j-1;
    end

    %Step N4
    if(mod(g(j)-z(j),2)==0)
        z(j)=z(j)+2;
        continue;
    end

    %Step N5
    t=z(j)-2;
    if(t<0)
        return
    end

    if(t<=z(j-1))
        t=t+2*(t<z(j-1))+1;
    end
    z(j)=t;
end

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
% batch input data setup:
isFullData = 1;
data_masterdir = GetCurrentDataDir();

const_ClusGroup = 2;
const_Cluster = 1; % This is all cells
M_fishset = GetFishStimset();
M_stimrange = GetStimRange();

range_fish =  1:18; % range_fish = GetFishRange();

%% custom params here:

%%
for i = 1:length(range_fish),
    i_fish = range_fish(i);
    disp(i_fish);
    fishset = M_fishset(i_fish);

    LoadFullFish(hfig,i_fish,isFullData);
    absIX = getappdata(hfig,'absIX');

    %% Cluster indexing
    i_ClusGroup = const_ClusGroup;% M_ClusGroup(i);
    i_Cluster = const_Cluster;% M_Cluster(i);
    cIX = LoadCluster_Direct(i_fish,i_ClusGroup,i_Cluster,absIX);

    %% get time index
    timelists = getappdata(hfig,'timelists');
%     timelists_names = getappdata(hfig,'timelists_names'); % use to screen
    stimrange = M_stimrange(i);

    tIX = GetTimeIndex_Direct(stimrange,timelists,fishset);

    %% get data matrix
    M = GetTimeIndexedData_Default_Direct(hfig,cIX,tIX);
    %M_0 = GetTimeIndexedData_Default_Direct(hfig,[],tIX,'isAllCells');

    %% ------custom code here---------

%         gIX = Kmeans_Direct(M,numK);

    %% (optional) save cluster
%     name = 'test';
%     clusgroupID = 2;
%     clusIDoverride = ; %Saved in Group2, Clusters 3&4
%     SaveCluster_Direct(cIX,gIX,absIX,i_fish,name,clusgroupID,clusIDoverride);

end

author: xiuyechen

repo: FishExplorer


text:
function [rsun,rmoon,gmst]=sunmoonpos(tutc,erpv)

tut=timeadd(tutc,erpv(3));

%calculate sun and moon position in ECI
[rsun,rmoon]=sunmoonpos_eci(tut);

%calculate the transition matrix from ECI to ECEF
[U,gmst]=eci2ecef(tutc,erpv);

%calculate sun and moon position in ECEF
rsun=U*rsun;
rmoon=U*rmoon;

return

repo: GINav


text:
function fx = p43_fun ( n, x )

%*****************************************************************************80
%
%% P43_FUN evaluates the integrand for problem 43.
%
%  Discussion:
%
%    The problem has a parameter ALPHA that can be set by calling
%    P43_PARAM_SET.
%
%    The suggested parameter range is 0.1 <= ALPHA <= 2.0.
%
%    The integrand has an algebraic endpoint singularity at X = 1
%    times a singular factor.
%
%  Interval:
%
%    0 <= x <= 1
%
%  Integrand:
%
%    ( ln ( 1 / x ) )^( alpha - 1 )
%
%  Exact Integral:
%
%    Gamma(alpha)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Robert Piessens, Elise de Doncker-Kapenga,
%    Christian Ueberhuber, David Kahaner,
%    QUADPACK: A Subroutine Package for Automatic Integration,
%    Springer, 1983, page 84.
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N), the evaluation points.
%
%    Output, real FX(N), the integrand values.
%
  alpha = p43_param_get ( );

  fx(1:n) = ( log ( 1.0 ./ x(1:n) ) ).^( alpha - 1.0 );

  i = find ( x == 0.0 | x == 1.0 );
  fx(i) = 0.0;

  return
end


text:
function [ fea, out ] = ex_linearelasticity4( varargin )
%EX_LINEARELASTICITY4 Stress calculation of an I-beam attached to two brackets.
%
%   [ FEA, OUT ] = EX_LINEARELASTICITY4( VARARGIN ) Example to calculate displacements and
%   stresses for an I-beam suppored by two brackets with circular holes.
%
%   Accepts the following property/value pairs.
%
%       Input       Value/{Default}        Description
%       -----------------------------------------------------------------------------------
%       E           scalar {200e9}         Modulus of elasticity
%       nu          scalar {0.3}           Poissons ratio
%       force       scalar {1e5}           Load force
%       l           scalar {0.4}           Length of I-beam
%       ilev        scalar {2}             Grid regfinement level
%       sfun        string {sflag1}        Shape function for displacements
%       iplot       scalar 0/{1}           Plot solution (=1)
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       fea         struct                 Problem definition struct
%       out         struct                 Output struct

% Copyright 2013-2022 Precise Simulation, Ltd.

cOptDef = { ...
  'E',        200e9; ...
  'nu',       0.3; ...
  'force',    1e5; ...
  'l',        0.4; ...
  'ilev',     2; ...
  'sfun',     'sflag1'; ...
  'iplot',    1; ...
  'tol',      0.42; ...
  'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid       = opt.fid;

% Geometry definition.
fea.sdim = { 'x' 'y' 'z' };   % Coordinate names.

% Grid generation.
fea.grid = get_grid( opt.ilev );

% Problem definition.
fea = addphys(fea,@linearelasticity);
fea.phys.el.eqn.coef{1,end} = { opt.nu };
fea.phys.el.eqn.coef{2,end} = { opt.E  };
fea.phys.el.sfun            = { opt.sfun opt.sfun opt.sfun };

% Boundary conditions.
dtol     = sqrt(eps);
fixbdr   = findbdr( fea, ['(sqrt(x.^2+z.^2)<=0.03+sqrt(eps))&(z>=-sqrt(eps))'] );
forcebdr = findbdr( fea, ['abs(y)>=0.2-sqrt(eps)'] );

% Fix boundaries (set zero Dirichlet BCs).
n_bdr  = max(fea.grid.b(3,:));        % Number of boundaries.
bctype = num2cell( zeros(3,n_bdr) );  % First set homogenous Neumann BCs everywhere.
[bctype{:,fixbdr}] = deal( 1 );       % Set Dirchlet BCs for right boundary.
fea.phys.el.bdr.coef{1,5} = bctype;

% Apply negative z-load to left boundary.
bccoef = num2cell( zeros(3,n_bdr) );
[bccoef{3,forcebdr}] = deal(-opt.force);
fea.phys.el.bdr.coef{1,end} = bccoef;

% Parse and solve problem.
fea       = parsephys( fea );
fea       = parseprob( fea );
fea.sol.u = solvestat( fea, 'fid', fid );

% Postprocessing.
if ( opt.iplot>0 )
  DSCALE = 5000;

  subplot(1,2,1)
  postplot( fea, 'surfexpr', 'sqrt(u^2+v^2+w^2)', 'linestyle', 'none' )
  title( 'Total displacement' )
  view([30 20])

  subplot(1,2,2)
  dp = zeros(size(fea.grid.p));
  for i=1:3
    dp(i,:) = DSCALE*evalexpr( fea.dvar{i}, fea.grid.p, fea );
  end
  fea_disp.grid   = fea.grid;
  fea_disp.grid.p = fea_disp.grid.p + dp;
  plotgrid( fea_disp )
  title(['Displacement plot (at ',num2str(DSCALE),' times scale)'])
  view([30 20])

end

% Error check.
disp_max_ref = 6.204e-6;
xdisp = fea.sol.u(fea.eqn.dofm{1}(:));
ydisp = fea.sol.u(fea.eqn.dofm{2}(:)+fea.eqn.ndof(1));
zdisp = fea.sol.u(fea.eqn.dofm{3}(:)+sum(fea.eqn.ndof(1:2)));
disp  = sqrt(xdisp.^2+ydisp.^2+zdisp.^2);
disp_max = max(disp);

svm_max_ref = 4.410e6;
svm = evalexpr( fea.phys.el.eqn.vars{1,2}, fea.grid.p, fea );
svm_max = max(svm);

out.disp_max = disp_max;
out.svm_max  = svm_max;
out.err(1)   = abs(disp_max - disp_max_ref)/abs(disp_max_ref);
out.err(2)   = abs(svm_max - svm_max_ref)/abs(svm_max_ref);
out.pass     = all(out.err<opt.tol);

if ( nargout==0 )
  clear fea out
end

%------------------------------------------------------------------------------%
function [ grid ] = get_grid( ilev )

n0 = 12;
nr = n0*2^(ilev-1);
r  = 0.03;   % Radius of bracket holes.
t  = 0.03;   % Thickness of brackets.
grid01 = ringgrid( 3*2^(ilev-1), 4*nr, r, r+t, [0;0] );
indc01 = selcells( grid01, 'y<=sqrt(eps)' );
grid01 = delcells( grid01, indc01 );

grid02 = holegrid( nr, 3*2^(ilev-1), (r+t)*[-1 1;-1 1], r, [0;0] );
indc02 = selcells( grid02, 'y>=-sqrt(eps)' );
grid02 = delcells( grid02, indc02 );
grid2d = gridmerge( grid01, findbdr(grid01,'y<=sqrt(eps)'), grid02, findbdr(grid02,'y>=-sqrt(eps)') );
t_br = 0.02;   % Width/depth of brackets.
d_br = 0.05;   % Separation distance between brackets.
grid1 = gridextrude( grid2d, 2^(ilev-1), t_br );
grid1 = gridrotate( grid1, pi/2, 1 );
grid2 = grid1;
grid1.p(2,:) = grid1.p(2,:) - d_br/2;
grid2.p(2,:) = grid2.p(2,:) + t_br + d_br/2;

% Create grids for the I-beam.
w_ib = 0.16;   % Beam width.
l_ib = 0.4;    % Beam length.
t_ib = 0.01;   % Beam thickness.
h_ib = 0.1;    % Beam height.

x_in    = linspace(-(r+t),r+t,n0+1);
x_coord = [ -w_ib/2 x_in w_ib/2];
y_coord = [ -0.2 -0.175 -0.15 -0.125 -0.1 -0.075 -(d_br/2+t_br) -d_br/2 0 d_br/2 d_br/2+t_br 0.075 0.1 0.125 0.15 0.175 0.2 ];
for i=2:ilev
  x_coord = sort([ x_coord [x_coord(1:end-1) + x_coord(2:end)]/2 ]);
  y_coord = sort([ y_coord [y_coord(1:end-1) + y_coord(2:end)]/2 ]);
end
grid3 = blockgrid( x_coord, y_coord, 2^(ilev-1), ...
                   [-w_ib/2 w_ib/2;-0.2 0.2;-(r+t)-t_ib -(r+t)] );
tm_ib = 2*(x_in(2)-x_in(1));
grid4 = blockgrid( 2*2^(ilev-1), y_coord, 5*2^(ilev-1), ...
                   [-tm_ib/2 tm_ib/2;-0.2 0.2;-(r+t)-t_ib-h_ib -(r+t)-t_ib] );
grid5 = grid3;
grid5.p(3,:) = grid5.p(3,:) - t_ib - h_ib;

% Merge grids.
tol = sqrt(eps)*1e3;
grid = gridmerge( grid1, findbdr(grid1,['z<=',num2str(-(r+t-tol))]), ...
                  grid3, findbdr(grid3,['z>=',num2str(-(r+t+tol))]) );
grid = gridmerge( grid2, findbdr(grid1,['z<=',num2str(-(r+t-tol))]), ...
                  grid,  findbdr(grid, ['(z>=',num2str(-(r+t+tol)),')&'...
                                        '(z<=',num2str(-(r+t-tol)),')']) );
grid = gridmerge( grid,  findbdr(grid,  ['z<=',num2str(-(r+t+t_ib-tol))]), ...
                  grid4, findbdr(grid4, ['z>=',num2str(-(r+t+t_ib+tol))]), 1 );
grid = gridmerge( grid,  findbdr(grid,  ['z<=',num2str(-(r+t+t_ib+h_ib-tol))]), ...
                  grid5, findbdr(grid5, ['z>=',num2str(-(r+t+t_ib+h_ib+tol))]), 2 );

author: precise-simulation

repo: featool-multiphysics


text:
function theta = initializeParameters(hiddenSize, visibleSize)

%% Initialize Parameters Randomly Based on Layer Sizes
r  = sqrt(6) / sqrt(hiddenSize+visibleSize+1);
W1 = rand(hiddenSize, visibleSize) * 2 * r - r;
W2 = rand(visibleSize, hiddenSize) * 2 * r - r;

b1 = zeros(hiddenSize, 1);
b2 = zeros(visibleSize, 1);

%% Compress Parameters into a Single Vector (theta)
theta = [W1(:) ; W2(:) ; b1(:) ; b2(:)];

end

author: zellyn


text:
function result  = gomp_mmv(Dict, X, K, L, T, epsilon, options)
    if(isobject(Dict))
        Dict = double(Dict);
    end
    if ~isa(Dict, 'double')
        error('Dict must be a double matrix');
    end
    if nargin < 4
        L = 1;
    end
    if nargin < 5
        T = 1;
    end
    if nargin < 6
        epsilon = 1e-3;
    end
    if nargin < 7
        options = struct;
    end
    if ~isfield(options, 'sparse_output')
        options.sparse_output = 1;
    end
    if ~isfield(options, 'verbose')
        options.verbose = 0;
    end
    result = mex_gomp_mmv(Dict, X, K, L, T, epsilon, ...
        options.sparse_output,...
        options.verbose);
end

author: indigits

repo: sparse-plex