text:
function [pvec, pstruct] = tapas_hgf_whatworld_transp(r, ptrans)
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2013 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Number of states whose contingencies have to be learned
ns = r.c_prc.n_states;

% Number of elements of the transition matrix
ntr = ns^2;

pvec    = NaN(1,length(ptrans));
pstruct = struct;

pvec(1:ntr)        = ptrans(1:ntr);                     % mu2_0
pstruct.mu2_0      = pvec(1:ntr);
pvec(ntr+1:2*ntr)  = exp(ptrans(ntr+1:2*ntr));          % sa2_0
pstruct.sa2_0      = pvec(ntr+1:2*ntr);
pvec(2*ntr+1)      = ptrans(2*ntr+1);                   % mu3_0
pstruct.mu3_0      = pvec(2*ntr+1);
pvec(2*ntr+2)      = exp(ptrans(2*ntr+2));              % sa3_0
pstruct.sa3_0      = pvec(2*ntr+2);
pvec(2*ntr+3)      = tapas_sgm(ptrans(2*ntr+3),r.c_prc.kaub); % ka
pstruct.ka         = pvec(2*ntr+3);
pvec(2*ntr+4)      = ptrans(2*ntr+4);                   % om
pstruct.om         = pvec(2*ntr+4);
pvec(2*ntr+5)      = tapas_sgm(ptrans(2*ntr+5),r.c_prc.thub); % th
pstruct.th         = pvec(2*ntr+5);

return;

author: translationalneuromodeling

repo: tapas


text:
function []=ps_load_initial_isce(data_inc)
%PS_LOAD_INITIAL Initial load of PS files into matlab workspace
%   PS_LOAD_INITIAL() loads various files and stores them
%   in various workspaces. The version number, PSVER, is set to 1.
%
%   Andy Hooper, June 2006
%
%   ======================================================================
%   06/2006 AH: Initialize day variable
%   07/2006 AH: Fix error in la and bperp matrix calculation
%   07/2006 AH: Add patch compatibility
%   09/2006 AH: Store wrapped phase in separate workspace
%   11/2007 AH: Change bperp to be the bperp at 0 m
%   01/2008 AH: Selection of bperp files made more rigourous
%   03/2009 AH: Change bperp back to be the bperp at local height
%   10/2009 AH: Add n_image to mat file
%   04/2010 KS: Fixed small issue with new logit function
%   10/2009 MA: Oversampling factor file introduced
%   12/2012 DB: Add compatiblility with Matlab2012B, keep backward compatible
%   01/2016 DB: Replace save with stamps_save which checks for var size when saving
%   04/2017 DB: Include fix for isce2stamps calamp.out slc filenaming
%   04/2018 DB: Add support for baseline grids
%   03/2022 DHTM: change master to reference and ingesting TomoSAR_interflist
%   ======================================================================

%NB IFGS assumed in ascending date order
fprintf('Loading data into matlab...\n')

phname=['pscands.1.ph'];            % for each PS candidate, a float complex value for each ifg
ijname=['pscands.1.ij'];            % ID# Azimuth# Range# 1 line per PS candidate
bperpname=['bperp.1.in'];           % in meters 1 line per ifg
dayname=['day.1.in'];               % YYYYMMDD, 1 line per ifg
masterdayname=['reference_day.1.in'];  % YYYYMMDD
llname=['pscands.1.ll'];            % 2 float values (lon and lat) per PS candidate
daname=['pscands.1.da'];            % 1 float value per PS candidate
hgtname=['pscands.1.hgt'];          % 1 float value per PS candidate
laname=['look_angle.1.in'];         % grid of look angle values
headingname=['heading.1.in'];       % satellite heading
lambdaname=['lambda.1.in'];         % wavelength
calname=['calamp.out'];             % amplitide calibrations
widthname=['width.txt'];            % width of interferograms
lenname=['len.txt'];                % length of interferograms

psver=1;
incsavename=['inc',num2str(psver)];
lasavename=['la',num2str(psver)];

% Get matlab version as function arguments change with the matlab version
matlab_version = version('-release');           % [DB] getting the matlab version
matlab_version = str2num(matlab_version(1:4));  % [DB] the year

if ~exist(dayname,'file')
    dayname= ['../',dayname];
end
day = load(dayname);

% ingesting the TomoSAR ComSAR name
TomoSAR_interflist = 'TomoSAR_interflist.mat';
if ~exist(TomoSAR_interflist,'file')
    TomoSAR_interflist= ['../',TomoSAR_interflist];
    if ~exist(TomoSAR_interflist,'file')
        % assume PSDS output is the same as PSI, no need to index-ing
        TomoSAR_ind = 1:size(day,1);
    else
        load(TomoSAR_interflist)
        [~,TomoSAR_ind]=ismember(TomoSAR_interflist,day);
    end
else
    load(TomoSAR_interflist)
    [~,TomoSAR_ind]=ismember(TomoSAR_interflist,day);
end

day = day(TomoSAR_ind);

year=floor(day/10000);
month=floor((day-year*10000)/100);
monthday=day-year*10000-month*100;
slave_day=datenum(year,month,monthday);
[slave_day,day_ix]=sort(slave_day);

if ~exist(masterdayname,'file')
    masterdayname= ['../',masterdayname];
end
master_day=load(masterdayname);
master_day_yyyymmdd=master_day;
year=floor(master_day/10000);
month=floor((master_day-year*10000)/100);
monthday=master_day-year*10000-month*100;
master_day=datenum(year,month,monthday);

master_ix=sum(slave_day<master_day)+1;
day=[slave_day(1:master_ix-1);master_day;slave_day(master_ix:end)]; % insert master day

%% bperp one value per slave
if ~exist(bperpname,'file')
    bperpname= ['../',bperpname];
end
bperp=load(bperpname);
bperp = bperp(TomoSAR_ind);

bperp=bperp(day_ix);
bperp=[bperp(1:master_ix-1);0;bperp(master_ix:end)]; % insert master-master bperp (zero)
n_ifg=size(bperp,1);
n_image=n_ifg;

%% heading
if ~exist(headingname,'file')
    headingname= ['../',headingname];
end
heading=load(headingname);
if isempty(heading)
    error('heading.1.in is empty')
end
setparm('heading',heading,1);

%% wavelength
if ~exist(lambdaname,'file')
    lambdaname= ['../',lambdaname];
end
lambda=load(lambdaname);
setparm('lambda',lambda,1);

%% Radar coordinates
ij=load(ijname);
n_ps=size(ij,1);

if ~exist(calname,'file')
    calname= ['../',calname];
end
if exist(calname,'file')
    [calfile,calconst]=textread(calname,'%s%f');
    caldate=zeros(length(calfile),1);
    for i = 1 : length(calfile)
        aa=strread(calfile{i},'%s','delimiter','/');

        if strcmpi(aa{end-1},'reference');
                bb=str2num(aa{end-2}(end-7:end));
        else
            bb=str2num(aa{end}(1:8));
            if isempty(bb)
                bb=str2num(aa{end-1}(1:8));
            end
        end
        caldate(i)=bb;
    end
    not_master_ix=caldate~=master_day_yyyymmdd;
    caldate=caldate(not_master_ix);
    calconst=calconst(not_master_ix);
    [Y,I]=sort(caldate);
    calconst=calconst(I);
else
    calconst=ones(n_ifg-1,1);
end

fid=fopen(phname,'r');
ph=zeros(n_ps,n_ifg-1,'single');
byte_count=n_ps*2;
for i=1:n_ifg-1
    [ph_bit,byte_count]=fread(fid,[(n_ps)*2,1],'float');
    ph_bit=single(ph_bit);
    ph(:,i)=complex(ph_bit(1:2:end),ph_bit(2:2:end));
end

ph=ph(:,day_ix);
zero_ph=sum(ph==0,2);
nonzero_ix=zero_ph<=1;       % if more than 1 phase is zero, drop node
ph=ph./repmat(calconst',n_ps,1);  % scale amplitudes
ph=[ph(:,1:master_ix-1),ones(n_ps,1),ph(:,master_ix:end)]; % insert zero phase master-master ifg

%% LON LAT
fid=fopen(llname,'r');
lonlat=fread(fid,[2,inf],'float');
lonlat=lonlat';
fclose(fid);

ll0=(max(lonlat)+min(lonlat))/2;
xy=llh2local(lonlat',ll0)*1000;
xy=xy';
sort_x=sortrows(xy,1);
sort_y=sortrows(xy,2);
n_pc=round(n_ps*0.001);
bl=mean(sort_x(1:n_pc,:)); % bottom left corner
tr=mean(sort_x(end-n_pc:end,:)); % top right corner
br=mean(sort_y(1:n_pc,:)); % bottom right  corner
tl=mean(sort_y(end-n_pc:end,:)); % top left corner

theta=(180-heading)*pi/180;
if theta>pi
    theta=theta-2*pi;
end

rotm=[cos(theta),sin(theta); -sin(theta),cos(theta)];
xy=xy';
xynew=rotm*xy; % rotate so that scene axes approx align with x=0 and y=0
if max(xynew(1,:))-min(xynew(1,:))<max(xy(1,:))-min(xy(1,:)) &...
   max(xynew(2,:))-min(xynew(2,:))<max(xy(2,:))-min(xy(2,:))
    xy=xynew; % check that rotation is an improvement
    disp(['Rotating by ',num2str(theta*180/pi),' degrees']);
end

xy=single(xy');
[dummy,sort_ix]=sortrows(xy,[2,1]); % sort in ascending y order
xy=xy(sort_ix,:);
xy=[[1:n_ps]',xy];
xy(:,2:3)=round(xy(:,2:3)*1000)/1000; % round to mm

ph=ph(sort_ix,:);
ij=ij(sort_ix,:);
ij(:,1)=1:n_ps;
lonlat=lonlat(sort_ix,:);

savename=['ps',num2str(psver)];
stamps_save(savename,ij,lonlat,xy,bperp,day,master_day,master_ix,n_ifg,n_image,n_ps,sort_ix,ll0,calconst,master_ix,day_ix);
save psver psver

phsavename=['ph',num2str(psver)];
% save(phsavename,'ph');
stamps_save(phsavename,ph);

if exist(daname,'file')
  D_A=load(daname);
  D_A=D_A(sort_ix);
  dasavename=['da',num2str(psver)];
  stamps_save(dasavename,D_A);
end

if exist(hgtname,'file')
    fid=fopen(hgtname,'r');
    hgt=fread(fid,[1,inf],'float');
    hgt=hgt';
    hgt=hgt(sort_ix);
    fclose(fid);
    hgtsavename=['hgt',num2str(psver)];
    stamps_save(hgtsavename,hgt);
end

if ~exist(widthname,'file')
    widthname= ['../',widthname];
end
width=load(widthname);

if ~exist(lenname,'file')
    lenname= ['../',lenname];
end
len=load(lenname);

% try to see if inc angle exist, fall back to la if needed
if nargin<1 | isempty(data_inc)
    clear data_inc
    % might be in dirs above when called from patch dir
    % inc angle
    if ~exist(incname,'file')
        incname= ['..',filesep,incname];
        if ~exist(incname,'file')
            incname= ['..',filesep,incname];
        end
    end

    if exist(incname,'file')
         data_inc = single(load_isce(incname));
    end
end
% checking if the inc angle exist if not try look angle
if exist('data_inc','var')
    % inc angle exist
    % getting the position of the PS
    IND = sub2ind(size(data_inc),ij(:,3)+1,ij(:,2)+1);
    % storing the data
    inc=data_inc(IND);
    inc = inc*pi./180;
    stamps_save(incsavename,inc)
else
    % trying look angle instead
     if ~exist(laname,'file')
        laname= ['..',filesep,laname];
        if ~exist(incname,'file')
            laname= ['..',filesep,laname];
            clear laname
        end
     end
     % laod the data if it exists
     if exist(laname,'file')
         data_la = single(load_isce(laname));

        % getting the position of the PS
        IND = sub2ind(size(data_la),ij(:,3)+1,ij(:,2)+1);

        % storing the data
        la=data_la(IND);
        la = la*pi./180;
        stamps_save(lasavename,la)
        clear data_la
     end
end

updir=0;
bperpdir=dir(['..' filesep 'baselineGRID_*']);
if isempty(bperpdir)
    bperpdir=dir(['..' filesep '..' filesep 'baselineGRID_*']);
    updir=1;
end
if length(bperpdir)>0
    % make a set of baseline excluding the master
    bperp_mat=zeros(n_ps,n_image-1,'single');
    counter =1;
    for i=setdiff([1:n_image],master_ix);
        bperp_fname=['..' filesep 'baselineGRID_' datestr(day(i),'yyyymmdd')];
        if updir==1
            bperp_fname=['..' filesep bperp_fname];
        end
        bperp_grid = load_isce(bperp_fname);

        % check if the baseline grid is a full IFG grid or smaller grid
        if size(bperp_grid)==[width, len]
            % getting the position of the PS
            IND = sub2ind(size(bperp_grid),ij(:,3)+1,ij(:,2)+1);
            % storing the data
            bp0_ps=bperp_grid(IND);
            if i==1
                fprintf('yes \n')
            end
        else
            if i==1
                fprintf('no \n')
                [gridX,gridY]=meshgrid(linspace(0,width,size(bperp_grid,1)),linspace(0,len,size(bperp_grid,2)));
            end
            bp0_ps=griddata_version_control(gridX,gridY,bperp_grid',ij(:,3),ij(:,2),'linear',matlab_version);               % [DB] fix matlab2012 version and older

        end
        bperp_mat(:,counter)=bp0_ps;
        counter = counter+1;
    end
else
    bperp_mat=repmat(single(bperp([1:master_ix-1,master_ix+1:end]))',n_ps,1);
end
bpsavename=['bp',num2str(psver)];
stamps_save(bpsavename,bperp_mat);

author: DinhHoTongMinh

repo: TomoSAR


text:
function [pointst,pointsf]=ridges(tfr,hat,t,method,trace);
%RIDGES	Extraction of ridges.
%	[POINTST,POINTSF]=RIDGES(TFR,HAT,T,METHOD,TRACE) extracts the
%	ridges of a time-frequency distribution. These ridges are some
%	particular sets of curves deduced from the stationary points of
%	their  reassignment operators.
%
%	TFR    : time-frequency representation
%	HAT    : complex matrix of the reassignment vectors.
%	T      : the time instant(s).
%	METHOD : the chosen representation (default: 'tfrrsp').
%	TRACE  : if nonzero, the progression of the algorithm is shown
%					   (default : 0).
%
%	POINTST,POINTSF are two vectors for the time and frequency
%	coordinates of the stationary points of the reassignment.
%	Therefore, PLOT(POINTST,POINTSF,'.') shows the squeleton of the
%	representation.
%
%	Example :
%	 sig=fmlin(128,0.1,0.4); g=tftb_window(21,'kaiser');
%	 h=tftb_window(47,'Kaiser'); t=1:2:127;
%	 figure(1), [tfr,rtfr,hat]=tfrrspwv(sig,t,128,g,h);
%	 ridges(tfr,hat,t,'tfrrspwv',1);
%	 figure(2), [tfr,rtfr,hat]=  tfrrsp(sig,t,128,h);
%	 ridges(tfr,hat,t,'tfrrsp',1);
%
%	See also : FRIEDMAN.

%	F. Auger, August 1994, December 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

if (nargin<2),
 error('at least 2 parameters required');
end;

[tfrrow,tfrcol]=size(tfr);
[hatrow,hatcol]=size(hat);

if (nargin==2),
 t=1:tfrcol; method='tfrrsp'; trace=0;
elseif (nargin==3),
 method='tfrrsp'; trace=0;
elseif (nargin==4),
 trace=0;
end;

[trow,tcol] = size(t);
if (trow~=1),
 error('T must only have one row');
elseif (tfrrow~=hatrow)|(tfrcol~=hatcol),
 error('TFR and HAT must have the same size');
end;

Nf=tfrrow; frequencies=(1:Nf)';
threshold=sum(sum(tfr))*0.5/(tfrrow*tfrcol);

pointst=[]; pointsf=[];

if trace, fprintf ('\nRidge extraction: \n'); end;

method=upper(method);
if strcmp(method,'TFRRPWV') | strcmp(method,'TFRRPMH'),
 for icol=1:tfrcol, ti=t(icol);
  if trace, disprog(icol,tfrcol,10); end;
  indices=find((tfr(:,icol)>threshold)&(hat(:,icol)-frequencies==0));
  nbindices=length(indices);
  if (nbindices>0),
   pointst=[pointst;ones(nbindices,1)*ti];
   pointsf=[pointsf;indices/(2.0*Nf)];
  end;
 end;
elseif strcmp(method,'TFRRSPWV'),
 for icol=1:tfrcol, ti=t(icol);
  if trace, disprog(icol,tfrcol,10); end;
  indices=find((real(hat(:,icol))-frequencies==0)&...
               (imag(hat(:,icol))-icol==0)&...
               (tfr(:,icol)>threshold));
  nbindices=length(indices);
  if (nbindices>0),
   pointst=[pointst;ones(nbindices,1)*ti];
   pointsf=[pointsf;indices/(2.0*Nf)];
  end;
 end;
elseif strcmp(method,'TFRRSP')|strcmp(method,'TYPE1')
 for icol=1:tfrcol, ti=t(icol);
  if trace, disprog(icol,tfrcol,10); end;
  indices=find((real(hat(:,icol))-frequencies==0)&...
               (imag(hat(:,icol))-icol==0)&...
               (tfr(:,icol)>threshold));
  nbindices=length(indices);
  if (nbindices>0),
   pointst=[pointst;ones(nbindices,1)*ti];
   pointsf=[pointsf;indices/Nf];
  end;
 end;
else
 error('unknown representation');
end;

if (nargout==0),
 clf;
 plot(pointst,pointsf,'.')
 axis([min(t) max(t) 0 0.5]);
end;

author: HeLiangHIT


text:
path('toolbox',path);
%% read pcd and show it
clear;clc;close all;
filename = '../data/2D/n.mat';
load(filename);
P.pts = M.verts;
clear('M');
npts = size(P.pts, 1);
[P.bbox, P.diameter] = GS.compute_bbox(P.pts);

%% call farthest_sampling_by_sphere
P.sample_radius = P.diameter*0.04;
[P.spls,P.corresp] = farthest_sampling_by_sphere(P.pts, P.sample_radius);

%% call connect_by_inherit_neigh
k=5;
kdtree = kdtree_build(P.pts);
P.neigh = zeros(npts, k);
for i = 1:npts
    P.neigh(i,:)  = kdtree_k_nearest_neighbors(kdtree, P.pts(i,:), k)';
end
kdtree_delete( kdtree );

P.spls_adj = connect_by_inherit_neigh(P.pts, P.spls, P.corresp, P.neigh);

%% call edge_collapse
options.collapse_order = 1;
[P.spls, P.skel_adj,P.corresp] = edge_collapse(P.pts, P.spls, P.corresp, P.neigh, P.spls_adj, options);

%% refine

%% draw
figure; movegui('northeast');set(gcf,'color','white');hold on;
plot3( P.spls(:,1), P.spls(:,2), P.spls(:,3), '.r', 'markersize', 5);
axis off; axis equal;set(gcf,'Renderer','OpenGL');
GS.plot_connectivity(P.spls, P.skel_adj, 2);
view3d zoom;

author: taiya

repo: cloudcontr

save_path: github-repos/MATLAB/taiya-cloudcontr


text:
function MimSaveListOfPointsAndValues(file_path, file_name, xc, yc, zc, value_list, coordinate_system, template_image)
    % MimSaveListOfPointsAndValues. Exports a list of coordinates and values
    %
    %     Syntax
    %     ------
    %
    %         MimSaveListOfPointsAndValues(file_path, file_name, xc, yc, zc, value_list, template_image, reporting)
    %
    %             file_path       is the path where the output file is to be
    %                             stored
    %             file_name       is the file name
    %             xc, yc, zx      are the coordinates in PTK cordinates
    %             coordinate_system  a MimCoordinateSystem enumeration
    %                             specifying the coordinate system to use
    %             template_image  A PTKImage providing voxel size and image size
    %                             parameters
    %             reporting       an object implementing CoreReportingInterface
    %                             for reporting progress and warnings
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    if nargin < 7
        reporting.Error('MimSaveListOfPointsAndValues:BadArguments', 'No coordinate_system parameter specified');
    end

    if ~isa(coordinate_system, 'MimCoordinateSystem')
        reporting.Error('MimSaveListOfPointsAndValues:BadArguments', 'coordinate_system parameter is not of type MimCoordinateSystem');
    end

    CoreDiskUtilities.CreateDirectoryIfNecessary(file_path);
    results_file_name = fullfile(file_path, file_name);
    file_handle = fopen(results_file_name, 'w');

    number_points = length(xc);

    for index = 1 : number_points
        dicom_coords = MimImageCoordinateUtilities.ConvertFromPTKCoordinates([xc(index), yc(index), zc(index)], coordinate_system, template_image);
        coord_x = dicom_coords(1);
        coord_y = dicom_coords(2);
        coord_z = dicom_coords(3);

        density = value_list(index);
        if any(density(:) == 0)
            disp('Warning: zero density');
        end

        output_string = sprintf('%6.6g,%6.6g,%6.6g,%6.6g\r\n', coord_x, coord_y, coord_z, density);
        fprintf(file_handle, regexprep(output_string, ' ', ''));
    end

    fclose(file_handle);
end

author: tomdoel

repo: pulmonarytoolkit


text:
function tc_meanAmpsPlusTcs(tc);
% tc_meanAmpsPlusTcs(tc);
%
% Provide a quick summary of ROI response
% to different conditions, by plotting both
% the mean time courses and mean amplitudes
% to each condition. Combines the
% routines from tc_plotMeanTrials and
% tc_barMeanAmplitudes.
%
%
% ras 03/05.

% delete existing large legend
otherAxes = findobj('Parent', gcf, 'Type', 'axes');
delete(otherAxes);

%%%%%%%%%%%%%%%%%%%
% mean amplitudes %
%%%%%%%%%%%%%%%%%%%
axes('Parent', tc.ui.plot, 'Position',[.1 .2 .35 .6]);
hold on

h1 = gca;

tc_barMeanAmplitudes(tc, h1, 'difference');  % Rory prefers it this way

%%%%%%%%%%%%%%%%%%%%%
% mean time courses %
%%%%%%%%%%%%%%%%%%%%%
h2 = axes('Parent', tc.ui.plot, 'Position', [.55 .2 .4 .6]);
tc_plotMeanTrials(tc, h2, 2);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [sol] = ls_solaat(rhs)

% Yin Zhang, 1997

global probData

% ----- check input rhs -----
if issparse(rhs)
   error('RHS vector must be in full matrix format.');
end;
m = length(probData.PERM);
if (max(size(rhs)) ~= m) || (min(size(rhs)) ~= 1)
   error('No symbolic factor or input sizes mismatch.');
end;

% ------ back solve ------
sol = zeros(m,1);
sol(probData.PERM) = ls_blkslv(probData.XLNZ,probData.XSUPER,probData.XLINDX,probData.LINDX,probData.LNZ0,rhs(probData.PERM));

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [rest] = mne_omit_first_line(str)
%
% [rest] = mne_omit_first_line(str)
%
% Omit the first line in a multi-line string (useful for handling
% error messages)
%

%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%
%   Revision 1.2  2006/04/23 15:29:40  msh
%   Added MGH to the copyright
%
%   Revision 1.1  2006/04/17 11:52:15  msh
%   Added coil definition stuff
%
%
me='MNE:mne_omit_first_line';

lf = findstr(10,str);
if isempty(lf)
    rest = str;
else
    rest = str(lf+1:size(str,2));
end

return;

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ssx, y2] = dofit(y, Qd)
%DOFIT Do constrained least squares fit and reduce data for subsequent fits

% Fit y to design matrix in null space
y2 = Qd' * y;            % rotate y into that space: predicted y value
ssx = norm(y2)^2;        % sum of squares explained by fit


text:
% function x = synsq_cwt_iw(Tx, fs, opt)
%
% Inverse Synchrosqueezing transform of Tx with associated
% frequencies in fs.  This implements Eq. 5 of [1].
%
% 1. E. Brevdo, N.S. Fučkar, G. Thakur, and H-T. Wu, "The
% Synchrosqueezing algorithm: a robust analysis tool for signals
% with time-varying spectrum," 2011.
%
% Input:
%   Tx, fs: See help synsq_cwt_fw
%   opt: options structure (see help synsq_cwt_fw)
%      opt.type: type of wavelet used in synsq_cwt_fw
%
%      other wavelet options (opt.mu, opt.s) should also match
%      those used in synsq_cwt_fw
%
% Output:
%   x: reconstructed signal
%
% Example:
%   [Tx,fs] = synsq_cwt_fw(t, x, 32); % Synchrosqueezing
%   Txf = synsq_filter_pass(Tx, fs, -Inf, 1); % Pass band filter
%   xf = synsq_cwt_iw(Txf, fs);  % Filtered signal reconstruction
%
%---------------------------------------------------------------------------------
%    Synchrosqueezing Toolbox
%    Authors: Eugene Brevdo (http://www.math.princeton.edu/~ebrevdo/)
%---------------------------------------------------------------------------------
function x = synsq_cwt_iw(Tx, fs, opt)
    if nargin<3, opt = struct(); end
    if ~isfield(opt, 'type'), opt.type = 'morlet'; end

    % Find the admissibility coefficient Cpsi
    Css = synsq_adm(opt.type, opt);

    [na, N] = size(Tx);

    % Integration
    % Due to linear discretization of integral in log(fs), this becomes
    % a simple normalized sum.
    x = real(1/Css * ones(1,na) * Tx);
end


text:
function varargout=lpcpf2cc(varargin)
%
% For calling details please see v_lpcpf2cc.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_lpcpf2cc.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: lpcpf2cc.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_lpcpf2cc(varargin{:});
else
    v_lpcpf2cc(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function varargout=wt(d,varargin)
%% Continous Wavelet Transform
% Creates a figure of wavelet power in units of
% normalized variance.
%
% USAGE: [wave,period,scale,coi,sig95]=wt(d[,params])
%
% d: a time series
% wave: the wavelet transform of d
% period: a vector of "Fourier" periods associated with wave
% scale: a vector of wavelet scales associated with wave
% coi: the cone of influence
%
% Settings: Pad: pad the time series with zeros?
% .         Dj: Octaves per scale (default: '1/12')
% .         S0: Minimum scale
% .         J1: Total number of scales
% .         Mother: Mother wavelet (default 'morlet')
% .         MaxScale: An easier way of specifying J1
% .         MakeFigure: Make a figure or simply return the output.
% .         BlackandWhite: Create black and white figures
% .         AR1: the ar1 coefficient of the series
% .              (default='auto' using a naive ar1 estimator. See ar1nv.m)
%
% Settings can also be specified using abbreviations. e.g. ms=MaxScale.
% For detailed help on some parameters type help wavelet.
%
%
% Example:
%      wt([0:200;sin(0:200)],'dj',1/20,'bw','maxscale',32)
%
% (C) Aslak Grinsted 2002-2014
%
% http://www.glaciology.net/wavelet-coherence

% -------------------------------------------------------------------------
%The MIT License (MIT)
%
%Copyright (c) 2014 Aslak Grinsted
%
%Permission is hereby granted, free of charge, to any person obtaining a copy
%of this software and associated documentation files (the "Software"), to deal
%in the Software without restriction, including without limitation the rights
%to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
%copies of the Software, and to permit persons to whom the Software is
%furnished to do so, subject to the following conditions:
%
%The above copyright notice and this permission notice shall be included in
%all copies or substantial portions of the Software.
%
%THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
%IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
%FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
%AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
%LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
%OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
%THE SOFTWARE.
%---------------------------------------------------------------------------

% ------validate and reformat timeseries.
[d,dt]=formatts(d);

n=size(d,1);
sigma2=var(d(:,2));

%----------default arguments for the wavelet transform-----------
Args=struct('Pad',1,...      % pad the time series with zeroes (recommended)
    'Dj',1/12, ...    % this will do 12 sub-octaves per octave
    'S0',2*dt,...    % this says start at a scale of 2 years
    'J1',[],...
    'Mother','Morlet', ...
    'MaxScale',[],...   %a more simple way to specify J1
    'MakeFigure',(nargout==0),...
    'AR1','auto');
Args=parseArgs(varargin,Args,{'BlackandWhite'});
if isempty(Args.J1)
    if isempty(Args.MaxScale)
        Args.MaxScale=(n*.17)*2*dt; %automaxscale
    end
    Args.J1=round(log2(Args.MaxScale/Args.S0)/Args.Dj);
end

if strcmpi(Args.AR1,'auto')
    Args.AR1=ar1nv(d(:,2));

    if any(isnan(Args.AR1))
        error('Automatic AR1 estimation failed. Specify it manually (use arcov or arburg).')
    end
end

%----------------::::::::---------- Analyze: ---------:::::::::::::-----------------

[wave,period,scale,coi] = wavelet(d(:,2),dt,Args.Pad,Args.Dj,Args.S0,Args.J1,Args.Mother);

t=d(:,1);
power = (abs(wave)).^2 ;        % compute wavelet power spectrum
signif = wave_signif(1.0,dt,scale,0,Args.AR1,-1,-1,Args.Mother);
sig95 = (signif')*(ones(1,n));  % expand signif --> (J+1)x(N) array
sig95 = power ./ (sigma2*sig95);
Yticks = 2.^(fix(log2(min(period))):fix(log2(max(period))));

if Args.MakeFigure
    H=imagesc(t,log2(period),log2(abs(power/sigma2)));%#ok,log2(levels));  %*** or use 'contourfill'
    %logpow=log2(abs(power/sigma2));
    %[c,H]=contourf(t,log2(period),logpow,[min(logpow(:)):.25:max(logpow(:))]);
    %set(H,'linestyle','none')

    clim=get(gca,'clim'); %center color limits around log2(1)=0
    clim=[-1 1]*max(clim(2),3);
    set(gca,'clim',clim)

    HCB=colorbar;
    set(HCB,'ytick',-7:7);
    barylbls=rats(2.^(get(HCB,'ytick')'));
    barylbls([1 end],:)=' ';
    barylbls(:,all(barylbls==' ',1))=[];
    set(HCB,'yticklabel',barylbls);

    set(gca,'YLim',log2([min(period),max(period)]), ...
        'YDir','reverse', ...
        'YTick',log2(Yticks(:)), ...
        'YTickLabel',num2str(Yticks'), ...
        'layer','top')
    %xlabel('Time')
    ylabel('Period')
    hold on

    [c,h] = contour(t,log2(period),sig95,[1 1],'k'); %#ok
    set(h,'linewidth',2)
    %plot(t,log2(coi),'k','linewidth',3)
    tt=[t([1 1])-dt*.5;t;t([end end])+dt*.5];
    hcoi=fill(tt,log2([period([end 1]) coi period([1 end])]),'w');
    set(hcoi,'alphadatamapping','direct','facealpha',.5)

    hold off
    set(gca,'box','on','layer','top');
end
varargout={wave,period,scale,coi,sig95};
varargout=varargout(1:nargout);

author: grinsted

repo: wavelet-coherence


text:
function result = pam(x,kclus,vtype,stdize,metric,silhplot)

%PAM is the Partitioning Around Medoids clustering algorithm.
% It returns a list representing a clustering of the data into kclus
% clusters based on the search for kclus representative objects or medoids among the observations of
% the data set.
%
% The algorithm is fully described in:
%   Kaufman, L. and Rousseeuw, P.J. (1990),
%   "Finding groups in data: An introduction to cluster analysis",
%   Wiley-Interscience: New York (Series in Applied Probability and
%   Statistics), ISBN 0-471-87876-6.
%
% Required input arguments:
%       x : Data matrix (rows = observations, columns = variables)
%           or Dissimilarity matrix (if number of columns equals 1)
%   kclus : The number of desired clusters
%   vtype : Variable type vector (length equals number of variables)
%           Possible values are 1  Asymmetric binary variable (0/1)
%                               2  Nominal variable (includes symmetric binary)
%                               3  Ordinal variable
%                               4  Interval variable
%          (if x is a dissimilarity matrix vtype is not required.)
%
% Optional input arguments:
%     stdize : standardise the variables given by the x-matrix
%              Possible values are 0 : no standardisation (default)
%                                  1 : standardisation by the mean
%                                  2 : standardisation by the median
%              (if x is a dissimilarity matrix, stdize is ignored)
%     metric : Metric to be used
%              Possible values are 'eucli' Euclidian (all interval variables, default)
%                                  'manha' Manhattan
%                                  'mixed' Mixed (not all interval variables, default)
%              (if x is a dissimilarity matrix, metric is ignored)
%   silhplot : draws picture
%              Possible values are 0 : do not create a silhouette plot (default)
%                                  1 : create a silhouette plot
%
% I/O:
%   result=pam(x,kclus,vtype,'eucli',silhplot)
%
% Example (subtracted from the referenced book)
%   load ruspini.mat
%   result=pam(ruspini,2,[4 4],1);
%
% The output of PAM is a structure containing:
%   result.dys        : dissimilarities (read row by row from the
%                       lower dissimilarity matrix)
%   result.metric     : metric used
%   result.number     : number of observations
%   result.ttd        : Average silhouette width per cluster
%   result.ttsyl      : Average silhouette width for dataset
%   result.idmed      : Id of medoid observations
%   result.obj        : Objective function at the first two iterations
%   result.ncluv      : Cluster membership for each observation
%   result.clusinf    : Matrix, each row gives numerical information for
%                       one cluster. These are the cardinality of the cluster
%                       (number of observations), the maximal and average
%                       dissimilarity between the observations in the cluster
%                       and the cluster's medoid, the diameter of the cluster
%                       (maximal dissimilarity between two observations of the
%                       cluster), and the separation of the cluster (minimal
%                       dissimilarity between an observation of the cluster
%                       and an observation of another cluster).
%   result.sylinf     : Matrix, with for each observation i the cluster to
%                       which i belongs, as well as the neighbor cluster of i
%                       (the cluster, not containing i, for which the average
%                       dissimilarity between its observations and i is minimal),
%                       and the silhouette width of i.
%   result.nisol      : Vector, with for each cluster specifying whether it is
%                       an isolated cluster (L- or L*-clusters) or not isolated.
%                       A cluster is an L*-cluster iff its diameter is smaller than
%                       its separation.  A cluster is an L-cluster iff for each
%                       observation i the maximal dissimilarity between i and any
%                       other observation of the cluster is smaller than the minimal
%                       dissimilarity between i and any observation of another cluster.
%                       Clearly each L*-cluster is also an L-cluster.
%
% And PAM will create the silhouette plot if silhplot equals 1 (an empty bar indicated by
%                       zero is a sparse between two clusters).
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at:
%              http://wis.kuleuven.be/stat/robust.html
%
% Written by Guy Brys (May 2006)

%Checking and filling in the inputs
res1=[];
if (nargin<2)
    error('Two input arguments required')
elseif ((nargin<3) & (size(x,2)~=1) & (size(x,1)~=1))
    error('Three input arguments required')
elseif (nargin<3)
    if (size(x,2)==1)
        x = x';
    end
    res1.metric = 'unknown';
    res1.disv = x;
    lookup=seekN(x);
    res1.number = lookup.numb;   %(1+sqrt(1+8*size(x,1)))/2;
    stdize = 0;
    silhplot = 0;
elseif (nargin<4)
    stdize = 0;
    silhplot = 0;
    if (sum(vtype)~=4*size(x,2))
        metric = 'mixed';
    else
        metric = 'eucli';
    end
elseif (nargin<5)
    silhplot = 0;
    if (sum(vtype)~=4*size(x,2))
        metric = 'mixed';
    else
        metric = 'eucli';
    end
elseif (nargin<6)
    silhplot = 0;
end

%Replacement of missing values
for i=1:size(x,1)
    A=find(isnan(x(i,:)));
    if (~(isempty(A)))
        for j=A
            valmisdat=0;
            for c=1:size(x,2)
                if (c~=j)
                    [a,b] = sort(x(:,c));
                    if (~isempty(b(find(a==x(i,c)))))
                        valmisdat=valmisdat+find(a==x(i,c));
                    end
                end
            end
            x(i,j)=prctile(x(isnan(x(:,j))==0,j),100*valmisdat/(size(x,1)*(size(x,2)-1)));
        end
    end
end

%Standardization
if ((stdize==1) & (strcmp(metric,'eucli')))
    x = ((x - repmat(mean(x),size(x,1),1))./(repmat(std(x),size(x,1),1)));
elseif ((stdize==2) & (strcmp(metric,'eucli')))
    x = ((x - repmat(median(x),size(x,1),1))./(repmat(mad(x),size(x,1),1)));
end

%Calculating the dissimilarities with daisy
if (isempty(res1))
    res1=daisy(x,vtype,metric);
end

%Actual calculations (the second for latter use with CLUSPLOT)
[dys,ttd,ttsyl,idmed,obj,ncluv,clusinf,sylinf,nisol]=pamc(res1.number,kclus,[0 res1.disv]');
dys=res1.disv(lowertouppertrinds(res1.number));

%Create a silhouetteplot
if (silhplot==1)
    Y=sylinf(:,3);
    Y1=flipdim(Y,1);
    whitebg([1 1 1]);
    % we calculate b="a but with a bar with length zero if the objects
    % are from another cluster"
    % and h="objects but with a 0 between 2 clusters"="g with a 0 if
    % it is a sparse between 2 clusters"
    a=flipdim(Y1,1);
    b=[];
    g=sylinf(:,4);
    f=sylinf(:,1)-1;
    for j=1:res1.number
        b(j+f(j))=a(j);
        h(j+f(j))=g(j);
    end
    b1=flipdim(b,2);
    h1=flipdim(h,2);
    % we use this b1 and h1 to plot the barh (instead of a and g)
    barh(b1,1);
    title 'Silhouette Plot of Pam' ;
    xlabel('Silhouette width');
    YT=1:res1.number+(sylinf(res1.number,1)-1);
    set(gca,'YTick',YT);
    set(gca,'YTickLabel',h1);
    axis([min([Y' 0]),max([Y' 0]),0.5,res1.number+0.5+f(res1.number)]);
elseif ((silhplot~=0) & (silhplot~=1) & (nargin==6))
    error('silhplot must equals 0 or 1')
end

%Putting things together
result = struct('dys',dys,'metric',res1.metric,'number',res1.number,...
    'ttd',ttd,'ttsyl',ttsyl,'idmed',idmed,'obj',obj,'ncluv',ncluv,...
    'clusinf',clusinf,'sylinf',sylinf,'nisol',nisol,'x',x);

%------------
%SUBFUNCTIONS

function dv = lowertouppertrinds(n)

dv=[];
for i=0:(n-2)
    dv = [dv cumsum(i:(n-2))+repmat(1+sum(0:i),1,n-i-1)];
end

%---
function outn = seekN(x)

ok=0;
numb=0;
k=size(x,2);
sums=cumsum(1:k);
for i=1:k
    if(sums(i)==k)

        numb=i+1;
        ok=1;
    end
end
outn=struct('numb',numb,'ok',ok);

author: marianux

repo: ecg-kit

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


text:
%% Machine Learning Online Class
%  Exercise 5 | Regularized Linear Regression and Bias-Variance
%
%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the
%  exercise. You will need to complete the following functions:
%
%     linearRegCostFunction.m
%     learningCurve.m
%     validationCurve.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%

%% Initialization
clear ; close all; clc

%% =========== Part 1: Loading and Visualizing Data =============
%  We start the exercise by first loading and visualizing the dataset.
%  The following code will load the dataset into your environment and plot
%  the data.
%

% Load Training Data
fprintf('Loading and Visualizing Data ...\n')

% Load from ex5data1:
% You will have X, y, Xval, yval, Xtest, ytest in your environment
load ('ex5data1.mat');

% m = Number of examples
m = size(X, 1);

% Plot training data
plot(X, y, 'rx', 'MarkerSize', 10, 'LineWidth', 1.5);
xlabel('Change in water level (x)');
ylabel('Water flowing out of the dam (y)');

fprintf('Program paused. Press enter to continue.\n');
pause;

%% =========== Part 2: Regularized Linear Regression Cost =============
%  You should now implement the cost function for regularized linear
%  regression.
%

theta = [1 ; 1];
J = linearRegCostFunction([ones(m, 1) X], y, theta, 1);

fprintf(['Cost at theta = [1 ; 1]: %f '...
         '\n(this value should be about 303.993192)\n'], J);

fprintf('Program paused. Press enter to continue.\n');
pause;

%% =========== Part 3: Regularized Linear Regression Gradient =============
%  You should now implement the gradient for regularized linear
%  regression.
%

theta = [1 ; 1];
[J, grad] = linearRegCostFunction([ones(m, 1) X], y, theta, 1);

fprintf(['Gradient at theta = [1 ; 1]:  [%f; %f] '...
         '\n(this value should be about [-15.303016; 598.250744])\n'], ...
         grad(1), grad(2));

fprintf('Program paused. Press enter to continue.\n');
pause;

%% =========== Part 4: Train Linear Regression =============
%  Once you have implemented the cost and gradient correctly, the
%  trainLinearReg function will use your cost function to train
%  regularized linear regression.
%
%  Write Up Note: The data is non-linear, so this will not give a great
%                 fit.
%

%  Train linear regression with lambda = 0
lambda = 0;
[theta] = trainLinearReg([ones(m, 1) X], y, lambda);

%  Plot fit over the data
plot(X, y, 'rx', 'MarkerSize', 10, 'LineWidth', 1.5);
xlabel('Change in water level (x)');
ylabel('Water flowing out of the dam (y)');
hold on;
plot(X, [ones(m, 1) X]*theta, '--', 'LineWidth', 2)
hold off;

fprintf('Program paused. Press enter to continue.\n');
pause;

%% =========== Part 5: Learning Curve for Linear Regression =============
%  Next, you should implement the learningCurve function.
%
%  Write Up Note: Since the model is underfitting the data, we expect to
%                 see a graph with "high bias" -- slide 8 in ML-advice.pdf
%

lambda = 0;
[error_train, error_val] = ...
    learningCurve([ones(m, 1) X], y, ...
                  [ones(size(Xval, 1), 1) Xval], yval, ...
                  lambda);

plot(1:m, error_train, 1:m, error_val);
title('Learning curve for linear regression')
legend('Train', 'Cross Validation')
xlabel('Number of training examples')
ylabel('Error')
axis([0 13 0 150])

fprintf('# Training Examples\tTrain Error\tCross Validation Error\n');
for i = 1:m
    fprintf('  \t%d\t\t%f\t%f\n', i, error_train(i), error_val(i));
end

fprintf('Program paused. Press enter to continue.\n');
pause;

%% =========== Part 6: Feature Mapping for Polynomial Regression =============
%  One solution to this is to use polynomial regression. You should now
%  complete polyFeatures to map each example into its powers
%

p = 8;

% Map X onto Polynomial Features and Normalize
X_poly = polyFeatures(X, p);
[X_poly, mu, sigma] = featureNormalize(X_poly);  % Normalize
X_poly = [ones(m, 1), X_poly];                   % Add Ones

% Map X_poly_test and normalize (using mu and sigma)
X_poly_test = polyFeatures(Xtest, p);
X_poly_test = bsxfun(@minus, X_poly_test, mu);
X_poly_test = bsxfun(@rdivide, X_poly_test, sigma);
X_poly_test = [ones(size(X_poly_test, 1), 1), X_poly_test];         % Add Ones

% Map X_poly_val and normalize (using mu and sigma)
X_poly_val = polyFeatures(Xval, p);
X_poly_val = bsxfun(@minus, X_poly_val, mu);
X_poly_val = bsxfun(@rdivide, X_poly_val, sigma);
X_poly_val = [ones(size(X_poly_val, 1), 1), X_poly_val];           % Add Ones

fprintf('Normalized Training Example 1:\n');
fprintf('  %f  \n', X_poly(1, :));

fprintf('\nProgram paused. Press enter to continue.\n');
pause;

%% =========== Part 7: Learning Curve for Polynomial Regression =============
%  Now, you will get to experiment with polynomial regression with multiple
%  values of lambda. The code below runs polynomial regression with
%  lambda = 0. You should try running the code with different values of
%  lambda to see how the fit and learning curve change.
%

lambda = 1;
[theta] = trainLinearReg(X_poly, y, lambda);

% Plot training data and fit
figure(1);
plot(X, y, 'rx', 'MarkerSize', 10, 'LineWidth', 1.5);
plotFit(min(X), max(X), mu, sigma, theta, p);
xlabel('Change in water level (x)');
ylabel('Water flowing out of the dam (y)');
title (sprintf('Polynomial Regression Fit (lambda = %f)', lambda));

figure(2);
[error_train, error_val] = ...
    learningCurve(X_poly, y, X_poly_val, yval, lambda);
plot(1:m, error_train, 1:m, error_val);

title(sprintf('Polynomial Regression Learning Curve (lambda = %f)', lambda));
xlabel('Number of training examples')
ylabel('Error')
axis([0 13 0 100])
legend('Train', 'Cross Validation')

fprintf('Polynomial Regression (lambda = %f)\n\n', lambda);
fprintf('# Training Examples\tTrain Error\tCross Validation Error\n');
for i = 1:m
    fprintf('  \t%d\t\t%f\t%f\n', i, error_train(i), error_val(i));
end

fprintf('Program paused. Press enter to continue.\n');
pause;

%% =========== Part 8: Validation for Selecting Lambda =============
%  You will now implement validationCurve to test various values of
%  lambda on a validation set. You will then use this to select the
%  "best" lambda value.
%

[lambda_vec, error_train, error_val] = ...
    validationCurve(X_poly, y, X_poly_val, yval);

close all;
plot(lambda_vec, error_train, lambda_vec, error_val);
legend('Train', 'Cross Validation');
xlabel('lambda');
ylabel('Error');

fprintf('lambda\t\tTrain Error\tValidation Error\n');
for i = 1:length(lambda_vec)
	fprintf(' %f\t%f\t%f\n', ...
            lambda_vec(i), error_train(i), error_val(i));
end

fprintf('Program paused. Press enter to continue.\n');
pause;

author: benoitvallon

repo: coursera-machine-learning


text:
function [Y, X, incrementalCertainty] = ...
    estimate_displacement_based_upon_vector_summation2d(solution,...
    filterDirection,...
    sizeImage,...
    numberOfFilterDirections,...
    dk,ck,...
    t11,t22,t12)
% ESTIMATE_DISPLACEMENT_BASED_UPON_VECTOR_SUMMATION2D Estimates a displacement field based upon vector summation
%
% [Y, X, incrementalCertainty] = ...
%    estimate_displacement_based_upon_vector_summation2d(solution,...
%                                                        filterDirection,...
%                                                        sizeImage,...
%                                                        numberOfFilterDirections,...
%                                                        dk,ck,...
%                                                        t11,t22,t12)
%
% INPUT ARGUMENTS
% solution                      - Solution to use (6-7)
% sizeImage                     - Image size
% dk                            - Phase-difference for k different filter
%                                 directions
% ck                            - Corresponding certainties for k different
%                                 filter directions
% filterDirection               - k filter directions
% numberOfFilterDirections      - Number of filter directions
% t11                           - Tensor element (1,1)
% t12                           - Tensor element (1,2)
% t22                           - Tensor element (2,2)
%
% OPTIONAL INPUT ARGUMENTS
% N/A
%
% OUTPUT ARGUMENTS
% Y	 							- Displacement along Y
% X 							- Displacement along X
% incrementalCertainty 			- Certainty of displacement

% Copyright (c) 2011 Daniel Forsberg
% danne.forsberg@outlook.com
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

X = zeros(sizeImage);
Y = zeros(sizeImage);
incrementalCertainty = zeros(sizeImage);

switch solution
    case 6 %************************************************
        for k = 1:numberOfFilterDirections
            ckprim(:,:,k) = ck(:,:,k).*((filterDirection{k}(1).*t11 + filterDirection{k}(2).*t12).*filterDirection{k}(1)...
                +(filterDirection{k}(1).*t12 + filterDirection{k}(2).*t22).*filterDirection{k}(2));
        end

        for k = 1:numberOfFilterDirections
            X = X+ckprim(:,:,k).*dk(:,:,k)*filterDirection{k}(2);
            Y = Y+ckprim(:,:,k).*dk(:,:,k)*filterDirection{k}(1);
            incrementalCertainty = incrementalCertainty+ckprim(:,:,k);
        end
    case 7 %************************************************
        for k = 1:numberOfFilterDirections
            X = X+ck(:,:,k).*dk(:,:,k)*filterDirection{k}(2);
            Y = Y+ck(:,:,k).*dk(:,:,k)*filterDirection{k}(1);
            incrementalCertainty = incrementalCertainty+ck(:,:,k);
        end
end

incrementalCertainty = incrementalCertainty + eps;
X = X./incrementalCertainty;
Y = Y./incrementalCertainty;

author: fordanic

repo: image-registration


text:
function nop(varargin)
    %NOP %NOP Do nothing
    %
    % NOP( ... )
    %
    % A do-nothing function for use as a placeholder when working with callbacks
    % or function handles.

    % Intentionally does nothing

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function err = nmf_eucl_dist(X,Y)

err = sum(sum((X-Y).^2));

author: andrewssobral

repo: lrslibrary


text:
function model = setFeedingFastingConstraints(model, feedingStatus,fastingValue,storageValue)
% This function sets constraints corresponding either to feeding (i.e.,
% storage reactions are set to have a lower bound of 0 while
% the upper bounds are set to >0 (storage is enabled) or to fasting (i.e.,
% storage reactions are set to have a lower bound of <0 (use of stored metabolites is enabled) while
% the upper bounds are set to 0 (storage is disabled).
%
% function model = setFeedingFastingConstraints(model, feedingStatus,fastingValue,storageValue)
%
% INPUT
% model             model structure
% feedingStatus     'feeding' = default bounds for diet uptake fastingValue for all diet exchanges, use diet script to refine, organ storage is turned on;
%                   'fasting' = all diet uptakes are closed, organ storage release is turned on
% fastingValue      default: fastingValue = -10;
% storageValue      default: storageValue = 10;
%
% OUTPUT
% model             with updated constraints
%
% Ines Thiele, 2015/2016
if ~exist('fastingValue','var')
    fastingValue = -10;
end
if ~exist('storageValue','var')
    storageValue = 10;
end

% load amino acids that can be stored
AAStorage;
storedInLiver={'nad(c)';'nadp(c)';'fad(c)';...
    'thmtp(c)';'pydam(c)';'pydx(c)';'pydxn(c)';...
    'coa(c)';'thf(c)';'btn(c)';'retinol(c)';'retfa(c)';...
    '11_cis_retfa(c)';'9_cis_retfa(c)';'25hvitd2(c)';'vitd3(c)';...
    'avite1(c)';'avite2(c)';'chol(c)';'fe3(c)';'phyQ(c)';...
    'glygn2(c)'};
StoreInKidney = {'ascb_L(c)';'chol(c)';};
StoreInMuscle ={'thmpp(c)';'25hvitd2(c)';'vitd3(c)';'chol(c)';'glygn2(c)'};

StoreInAdi={'25hvitd2(c)';'vitd3(c)';'avite1(c)';'avite2(c)';'phyQ(c)';...
    'hdca(c)';'tmndnc(c)';'lnlc(c)';'tag_hs(c)';...
    'c226coa(c)';'doco13ecoa(c)';'lnlccoa(c)';'lnlncacoa(c)';'lnlncgcoa(c)';...
    'odecoa(c)';'pmtcoa(c)';'stcoa(c)';'tmndnccoa(c)'
    };
if isfield(model,'rxnGeneMat')
    model = rmfield(model,'rxnGeneMat');
end
% set all sinks to 0
for i = 1 : length(model.rxns)
    if length(strfind(model.rxns{i},'sink_'))==1 ...
            && length(strfind(model.rxns{i},'sink_pre_prot(r)'))==0  ...
            && length(strfind(model.rxns{i},'sink_Ser_Gly_Ala_X_Gly(r)'))==0 ...
            && length(strfind(model.rxns{i},'sink_5hpet(c)'))==0
        %&& length(strfind(model.rxns{i},'sink_Tyr_ggn(c)'))==0 % ...
        %% && length(strfind(model.rxns{i},'sink_citr(c)'))==0
        % && length(strfind(model.rxns{i},'sink_Ser_Gly_Ala_X_Gly(r)'))==0 ...
        if isfield(model,'Microbiota') && model.Microbiota(i) ==0 %no microbe sink
            model.lb(i)=0;
        elseif ~isfield(model,'Microbiota')
            model.lb(i)=0;
        end
    end
end

DMs= (find(~cellfun(@isempty,strfind(model.rxns,'DM_'))));
model.lb(DMs) = 0;

% but not at the beginning of the abbr --> leaves in the
% microbe sinks
if 1
    tmp = strmatch('sink_',model.rxns);
    model.lb(tmp)=-10;
end

for i = 1 : length(model.rxns)
    if strfind(model.rxns{i},'sink_')
        model.rxns{i} = regexprep(model.rxns{i},'\[c\]','(c)');
        model.rxns{i} = regexprep(model.rxns{i},'\[r\]','(r)');
        %    model.rxns{i}
    end
end
modelexchanges1 = strmatch('Diet_EX_',model.rxns);
modelexchanges2 = strmatch('Diet_Ex_',model.rxns);
modelexchanges = [modelexchanges1;modelexchanges2];

if strcmp(feedingStatus,'feeding') % storage no sinks
    model.lb(modelexchanges)=fastingValue;
    model.ub(modelexchanges)=0;
elseif strcmp(feedingStatus,'fasting') % storage no sinks
    model.lb(modelexchanges)=0;
    model.ub(modelexchanges)=0;
    % exception for water

end

if strcmp(feedingStatus,'feeding') % storage no sinks
    for i = 1 : length(storedInLiver)
        storedInLiver{i};
        rxnName = strcat('Liver_sink_',storedInLiver{i});
        %   rxnName
        model = changeRxnBounds(model,rxnName,0,'l');
        model = changeRxnBounds(model,rxnName, storageValue,'u');
    end
    for i = 1 : length(StoreInKidney)
        StoreInKidney{i};
        rxnName = strcat('Kidney_sink_',StoreInKidney{i});
        %  rxnName
        model = changeRxnBounds(model,rxnName,0,'l');
        model = changeRxnBounds(model,rxnName, storageValue,'u');
    end
    for i = 1 : length(StoreInMuscle)
        rxnName = strcat('Muscle_sink_',StoreInMuscle{i});
        %  rxnName
        model = changeRxnBounds(model,rxnName,0,'l');
        model = changeRxnBounds(model,rxnName, storageValue,'u');
    end
    for i = 1 : length(storageAA)
        rxnName = strcat('Muscle_',storageAA{i});
        %    rxnName
        model = changeRxnBounds(model,rxnName,0,'l');
        model = changeRxnBounds(model,rxnName, storageValue,'u');
    end
    for i = 1 : length(StoreInAdi)
        rxnName = strcat('Adipocytes_sink_',StoreInAdi{i});
        %   rxnName
        StoreInAdi{i};
        model = changeRxnBounds(model,rxnName,0,'l');
        model = changeRxnBounds(model,rxnName, storageValue,'u');
    end
    rxnName = 'Retina_sink_crvnc(c)';
    model = changeRxnBounds(model,rxnName,0,'l');
    model = changeRxnBounds(model,rxnName, storageValue,'u');
    rxnName = 'Heart_sink_chol(c)';
    model = changeRxnBounds(model,rxnName,0,'l');
    model = changeRxnBounds(model,rxnName, storageValue,'u');
    rxnName = 'Brain_sink_crvnc(c)';
    model = changeRxnBounds(model,rxnName,0,'l');
    model = changeRxnBounds(model,rxnName, storageValue,'u');
    rxnName = 'Brain_sink_chol(c)';
    model = changeRxnBounds(model,rxnName,0,'l');
    model = changeRxnBounds(model,rxnName, storageValue,'u');
    rxnName = 'RBC_sink_glygn2(c)';
    model = changeRxnBounds(model,rxnName,0,'l');
    model = changeRxnBounds(model,rxnName, storageValue,'u');
    rxnName = 'Skin_sink_vitd3(c)';
    model = changeRxnBounds(model,rxnName,0,'l');
    model = changeRxnBounds(model,rxnName, storageValue,'u');
elseif strcmp(feedingStatus,'fasting')
    for i = 1 : length(storedInLiver)
        rxnName = strcat('Liver_sink_',storedInLiver{i});
        model = changeRxnBounds(model,rxnName,fastingValue,'l');
        model = changeRxnBounds(model,rxnName,0,'u');
    end
    for i = 1 : length(StoreInKidney)
        rxnName = strcat('Kidney_sink_',StoreInKidney{i});
        model = changeRxnBounds(model,rxnName,fastingValue,'l');
        model = changeRxnBounds(model,rxnName,0,'u');
    end
    for i = 1 : length(StoreInMuscle)
        rxnName = strcat('Muscle_sink_',StoreInMuscle{i});
        model = changeRxnBounds(model,rxnName,fastingValue,'l');
        model = changeRxnBounds(model,rxnName,0,'u');
    end
    for i = 1 : length(storageAA)
        rxnName = strcat('Muscle_',storageAA{i});
        model = changeRxnBounds(model,rxnName,fastingValue,'l');
        model = changeRxnBounds(model,rxnName,0,'u');
    end
    for i = 1 : length(StoreInAdi)
        rxnName = strcat('Adipocytes_sink_',StoreInAdi{i});
        model = changeRxnBounds(model,rxnName,fastingValue,'l');
        model = changeRxnBounds(model,rxnName,0,'u');
    end
    rxnName = 'Retina_sink_crvnc(c)';
    model = changeRxnBounds(model,rxnName,fastingValue,'l');
    model = changeRxnBounds(model,rxnName,0,'u');
    rxnName = 'Heart_sink_chol(c)';
    model = changeRxnBounds(model,rxnName,fastingValue,'l');
    model = changeRxnBounds(model,rxnName,0,'u');
    rxnName = 'Brain_sink_crvnc(c)';
    model = changeRxnBounds(model,rxnName,fastingValue,'l');
    model = changeRxnBounds(model,rxnName,0,'u');
    rxnName = 'Brain_sink_chol(c)';
    model = changeRxnBounds(model,rxnName,fastingValue,'l');
    model = changeRxnBounds(model,rxnName,0,'u');
    rxnName = 'RBC_sink_glygn2(c)';
    model = changeRxnBounds(model,rxnName,fastingValue,'l');
    model = changeRxnBounds(model,rxnName,0,'u');
    rxnName = 'Skin_sink_vitd3(c)';
    model = changeRxnBounds(model,rxnName,fastingValue,'l');
    model = changeRxnBounds(model,rxnName,0,'u');
end

model.SetupInfo.FeedingStatus = feedingStatus;

author: opencobra

repo: cobratoolbox


text:
%% Mesh Smoothing and Optimization
%
%% Improve geometric mesh quality
%
% The function [node,elem] = optmesh(node,elem) will optimize the shape
% regularity of triangles in the input mesh (node,elem) and outputs a
% better mesh (node,elem).

clear all; close all;
load airfoilperturb
figure(1); subplot(1,2,1);
showmesh(node,elem); title('original mesh');
figure(2); subplot(1,2,1);
showmeshquality(node,elem); axis([0 1 0 2700]);
[node,elem] = optmesh(node,elem);
figure(1); subplot(1,2,2);
showmesh(node,elem); title('smoothed mesh');
figure(2); subplot(1,2,2);
showmeshquality(node,elem); axis([0 1 0 2700]);
%%
% We explain algorithms implemented in optimesh.m in the following.

%% ODT-based mesh smoothing
%
% In the function <../../../mesh/html/meshsmoothing.html meshsmoothing>, we move one node at a time inside its
% patch, which consists of all triangles surrounding this node, such that
% the interpolation error to a quadratic function is minimized. The
% function |meshsmoothing| will keep the topology of the input mesh, i.e.,
% the node index and connectivity of nodes are unchanged.
%
% In the simplest case, the scheme is to move the node to the average of
% circumenters of triangles in the local patch. Details can be found in the
% paper <http://math.uci.edu/~chenlong/CH2008.html ODTmesh>.
theta = [-2*pi/3 -pi/3 0 pi/3 2*pi/3 pi]';
node = [cos(theta), sin(theta)];
node(end+1,:) = 0;
elem = delaunayn(node);
node(end,:) = rand(1,2)*0.4;
figure(1); subplot(1,2,1);
showmesh(node,elem); findnode(node,'all','noindex');
c = circumcenter(node,elem);
hold on; plot(c(:,1),c(:,2),'r.','MarkerSize',16)
node(end,:) = mean(c);
plot(node(end,1),node(end,2),'b.','MarkerSize',16)
figure(1); subplot(1,2,2);
showmesh(node,elem); findnode(node,'all','noindex');

repo: Aero-matlab


text:
function sim = gaussianKernel(x1, x2, sigma)
%RBFKERNEL returns a radial basis function kernel between x1 and x2
%   sim = gaussianKernel(x1, x2) returns a gaussian kernel between x1 and x2
%   and returns the value in sim

% Ensure that x1 and x2 are column vectors
x1 = x1(:); x2 = x2(:);

% You need to return the following variables correctly.
sim = 0;

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return the similarity between x1
%               and x2 computed using a Gaussian kernel with bandwidth
%               sigma
%
%

sim = exp(-sum((x1-x2).^2)/(sigma.^2*2));

% =============================================================

end

author: Ayatans

repo: Machine-Learning-homework


text:
function x = p21_start ( n )

%*****************************************************************************80
%
%% P21_START returns a starting point for optimization for problem 21.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 March 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables X.
%
%    Output, real X(N), a starting point for the optimization.
%
  x = ones ( n, 1 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ttm] = tt_Fd_mtx2(tt_a, bound1, bound2, eps)
%TT-representation of the diffusion matrix
%   [TTM] = TT_FD_MTX2(TT_A, BOUND1, BOUND2, EPS) Computes TT
%   representation of a simplest discretization of the diffusion operator
%   with operator given in the QTT-format (TT_A).
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

d = size(tt_a, 1);
n = zeros(1,d);
for q=1:d
    n(q)=size(tt_a{q}, 1);
end;

if (max(size(bound1))==1) bound1 = bound1*ones(1,d); end;
if (max(size(bound2))==1) bound2 = bound2*ones(1,d); end;

ranks_a = tt_ranks(tt_a);

for q=1:d
    cur_tt = cell(d,1);

    cur_tt{1}=zeros(n(1),n(1),ranks_a(1));
    for k=1:ranks_a(1)
            lp1 = diag(-1*ones(n(1)-1,1), [-1])+diag(2*ones(n(1), 1), [0])+diag(-1*ones(n(1)-1,1), [1]);
            Mp1 = diag((1/6)*ones(n(1)-1,1), [-1])+diag((4/6)*ones(n(1), 1), [0])+diag((1/6)*ones(n(1)-1,1), [1]);
            if (bound1(1)==1) lp1(1,1)=1; end;
            if (bound2(1)==1) lp1(n(1),n(1))=1; end;
            h = 1/(n(1)+1);
            if ((bound1(1)==1)||(bound2(1)==1)) h = 1/n(1); end;
            if ((bound1(1)==1)&&(bound2(1)==1)) h = 1/(n(1)-1); end;
            Mlp1 = diag(tt_a{1}(2:n(1),k), [-1])+diag(tt_a{1}(2:n(1),k), [1]);
            Mlp1 = Mlp1 + diag(tt_a{1}(1:n(1),k)+[tt_a{1}(2:n(1),k)' tt_a{1}(n(1),k)]', [0])*0.5;
        if (q==1)
            cur_tt{1}(:,:,k) = Mlp1.*(lp1/h^2);
        else
            cur_tt{1}(:,:,k) = Mlp1.*Mp1;
%             cur_tt{1}(:,:,k) = diag(tt_a{1}(1:n(1),k)+[tt_a{1}(2:n(1),k)' tt_a{1}(n(1),k)]', [0])*0.5;
        end;
    end;

    for p=2:d-1
        cur_tt{p}=zeros(n(p),n(p),ranks_a(p-1),ranks_a(p));
        for k1=1:ranks_a(p-1)
            for k2=1:ranks_a(p)
                    lp1 = diag(-1*ones(n(p)-1,1), [-1])+diag(2*ones(n(p), 1), [0])+diag(-1*ones(n(p)-1,1), [1]);
                    Mp1 = diag((1/6)*ones(n(1)-1,1), [-1])+diag((4/6)*ones(n(1), 1), [0])+diag((1/6)*ones(n(1)-1,1), [1]);
                    if (bound1(p)==1) lp1(1,1)=1; end;
                    if (bound2(p)==1) lp1(n(p),n(p))=1; end;
                    h = 1/(n(p)+1);
                    if ((bound1(p)==1)||(bound2(p)==1)) h = 1/n(p); end;
                    if ((bound1(p)==1)&&(bound2(p)==1)) h = 1/(n(p)-1); end;
                    Mlp1 = diag(tt_a{p}(2:n(p),k1,k2), [-1])+diag(tt_a{p}(2:n(p),k1,k2), [1]);
                    Mlp1 = Mlp1 + diag(tt_a{p}(1:n(p),k1,k2)+[tt_a{p}(2:n(p),k1,k2)' tt_a{p}(n(p),k1,k2)]', [0])*0.5;
                if (p==q)
                    cur_tt{p}(:,:,k1,k2) = Mlp1.*(lp1/h^2);
                else
                    cur_tt{p}(:,:,k1,k2) = Mlp1.*Mp1;
%                     cur_tt{p}(:,:,k1,k2) = diag(tt_a{p}(1:n(p),k1,k2)+[tt_a{p}(2:n(p),k1,k2)' tt_a{p}(n(p),k1,k2)]', [0])*0.5;
                end;
            end;
        end;
    end;

    cur_tt{d}=zeros(n(d),n(d),ranks_a(d-1));
    for k=1:ranks_a(d-1)
            lp1 = diag(-1*ones(n(d)-1,1), [-1])+diag(2*ones(n(d), 1), [0])+diag(-1*ones(n(d)-1,1), [1]);
            Mp1 = diag((1/6)*ones(n(1)-1,1), [-1])+diag((4/6)*ones(n(1), 1), [0])+diag((1/6)*ones(n(1)-1,1), [1]);
            if (bound1(d)==1) lp1(1,1)=1; end;
            if (bound2(d)==1) lp1(n(d),n(d))=1; end;
            h = 1/(n(d)+1);
            if ((bound1(d)==1)||(bound2(d)==1)) h = 1/n(d); end;
            if ((bound1(d)==1)&&(bound2(d)==1)) h = 1/(n(d)-1); end;
            Mlp1 = diag(tt_a{d}(2:n(d),k), [-1])+diag(tt_a{d}(2:n(d),k), [1]);
            Mlp1 = Mlp1 + diag(tt_a{d}(1:n(d),k)+[tt_a{d}(2:n(d),k)' tt_a{d}(n(d),k)]', [0])*0.5;
        if (q==d)
            cur_tt{d}(:,:,k) = Mlp1.*(lp1/h^2);
        else
            cur_tt{d}(:,:,k) = Mlp1.*Mp1;
%             cur_tt{d}(:,:,k) = diag(tt_a{d}(1:n(d),k)+[tt_a{d}(2:n(d),k)' tt_a{d}(n(d),k)]', [0])*0.5;
        end;
    end;

    if (q==1)
        ttm = cur_tt;
    else
        ttm = ttm_add(ttm, cur_tt);
        ttm = tt_mat_compr(ttm, eps);
    end;
end;

end

author: oseledets

repo: TT-Toolbox


text:
function x = tuple_next_fast ( m, n, rank )

%*****************************************************************************80
%
%% TUPLE_NEXT_FAST computes the next element of a tuple space, "fast".
%
%  Discussion:
%
%    The elements are N vectors.  Each entry is constrained to lie
%    between 1 and M.  The elements are produced one at a time.
%    The first element is
%      (1,1,...,1)
%    and the last element is
%      (M,M,...,M)
%    Intermediate elements are produced in lexicographic order.
%
%    This code was written as a possibly faster version of TUPLE_NEXT.
%
%  Example:
%
%    N = 2,
%    M = 3
%
%    INPUT        OUTPUT
%    -------      -------
%    Rank          X
%    ----          ----
%   -1            -1 -1
%
%    0             1  1
%    1             1  2
%    2             1  3
%    3             2  1
%    4             2  2
%    5             2  3
%    6             3  1
%    7             3  2
%    8             3  3
%    9             1  1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 August 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the maximum entry in each component.
%    M must be greater than 0.
%
%    Input, integer N, the number of components.
%    N must be greater than 0.
%
%    Input, integer RANK, indicates the rank of the tuples.
%    Typically, 0 <= RANK < N**M; values greater than this are
%    legal and meaningful, being equivalent to the corresponding
%    value mod N**M.  RANK < 0 indicates that this is the first
%    call for the given values of (M,N).  Initialization is done,
%    and X is set to a dummy value.
%
%    Output, integer X(N), the next tuple of the given rank,
%    or a dummy value if initialization is being done.
%
  global tuple_next_fast_BASE

  if ( rank < 0 )

    if ( m <= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'TUPLE_NEXT_FAST - Fatal error!\n' );
      fprintf ( 1, '  M <= 0 is illegal.\n' );
      fprintf ( 1, '  M = %d\n', m );
      error ( 'TUPLE_NEXT_FAST - Fatal error!' );
    end

    if ( n <= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'TUPLE_NEXT_FAST - Fatal error!\n' );
      fprintf ( 1, '  N <= 0 is illegal.\n' );
      fprintf ( 1, '  N = %d\n', n );
      error ( 'TUPLE_NEXT_FAST - Fatal error!' );
    end

    tuple_next_fast_BASE(n) = 1;
    for i = n-1 : -1 : 1
      tuple_next_fast_BASE(i) = tuple_next_fast_BASE(i+1) * m;
    end

    x(1:n) = -1;

  else

    x(1:n) = mod ( floor ( rank ./ tuple_next_fast_BASE(1:n) ), m ) + 1;

  end

  return
end


text:
function B = fast_unitary_tranform(A, D, type)
% B = fast_unitary_tranform(A, D, type)
%
% Applies a fast unitary transform (DCT, DHT and WHT) to D * A
% after padding with zeros.
%
% type is either 'DCT', 'DHT' or 'WHT'.
%
% 6-December 2009, Version 1.3
% Copyright (C) 2009, Haim Avron and Sivan Toledo.

author: Sable

repo: mcbench-benchmarks


text:
function value = year_length_months_common ( y )

%*****************************************************************************80
%
%% YEAR_LENGTH_MONTHS_COMMON returns the number of months in a Common year.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year to be checked.
%
%    Output, integer VALUE, the number of months in the year.
%
  value = 12;

  return
end


text:
%%********************************************************
%%
%% [A,B] = qapread('~/ZXY/PAL/QAP/qapdata/nug12.dat');
%%
%% SDPNAL:
%% Copyright (c) 2008 by
%% Xinyuan Zhao, Defeng Sun, and Kim-Chuan Toh
%%********************************************************

  function [A,B] = qapread(fname);

  fid = fopen(fname,'r');
  if (fid == -1); error('file cannot be opened'); end

  [datavec,count] = fscanf(fid,'%c');
  fclose('all');
  linefeeds = findstr(datavec,char(10));
  datavec(linefeeds) = blanks(length(linefeeds));
  datavec = sscanf(datavec,'%f');

  n = datavec(1); n2 = n*n;
  aa = datavec(2:n2+1);
  bb = datavec(n2+2:2*n2+1);

  A = reshape(aa,n,n);
  B = reshape(bb,n,n);
%%********************************************************

author: intellhave

repo: SDRSAC

save_path: github-repos/MATLAB/intellhave-SDRSAC


text:
function a = ref_random ( m, n, prob, key )

%*****************************************************************************80
%
%% REF_RANDOM returns a random row echelon matrix.
%
%  Definition:
%
%    1) the first nonzero entry in any row is 1.
%
%    2) the first nonzero entry in row I occurs in a later column
%       than the first nonzero entry of every previous row.
%
%    3) rows that are entirely zero occur after all rows with
%       nonzero entries.
%
%  Example:
%
%    M = 6, N = 5, PROB = 0.8
%
%     1.0  0.3  0.2  0.0  0.5
%     0.0  0.0  1.0  0.7  0.9
%     0.0  0.0  0.0  1.0  0.3
%     0.0  0.0  0.0  0.0  1.0
%     0.0  0.0  0.0  0.0  0.0
%     0.0  0.0  0.0  0.0  0.0
%
%  Properties:
%
%    A is generally not symmetric: A' /= A.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the order of the matrix.
%
%    Input, real PROB, the probability that the 1 in the next
%    row will be placed as early as possibly.
%    Setting PROB = 1 forces the 1 to occur immediately, setting
%    PROB = 0 forces the entire matrix to be zero.  A more reasonable
%    value might be PROB = 0.8 or 0.9.
%
%    Input, integer KEY, a positive value that selects the data..
%
%    Output, real A(M,N), the matrix.
%
  a = zeros ( m, n );

  jprev = 0;

  seed = key;

  for i = 1 : m

    jnew = 0;

    for j = 1 : n

      if ( j <= jprev )
        a(i,j) = 0.0;
      elseif ( jnew == 0 )
        [ temp, seed ] = r8_uniform_01 ( seed );
        if ( temp <= prob )
          jnew = j;
          a(i,j) = 1.0;
        else
          a(i,j) = 0.0;
        end
      else
        [ a(i,j), seed ] = r8_uniform_01 ( seed );
      end

    end

    if ( jnew == 0 )
      jnew = n + 1;
    end

    jprev = jnew;

  end

  return
end


text:
function UFexport (list, check, tmp)
%UFEXPORT export to Matrix Market and Rutherford/Boeing formats
%
% Example:
%   UFexport ;          % export the entire collection
%   UFexport (list) ;   % just export matrices whose id's are given in the list
%   UFexport (list, 'check') ;      % also read them back in, to check
%
% If the list is empty, all matrices in the collection are exported.
% A 3rd argument tmp changes the tmp directory for UFread.
%
% See also UFget, UFwrite, RBio, mwrite.

% Copyright 2006-2007, Timothy A. Davis

%-------------------------------------------------------------------------------
% get the input arguments
%-------------------------------------------------------------------------------

index = UFget ;
nmat = length (index.nrows) ;

if (nargin < 1 || isempty (list))
    list = 1:nmat ;
end

check = ((nargin > 1) && strcmp (check, 'check')) ;

if (nargin < 3)
    tmp = '' ;
end

%-------------------------------------------------------------------------------
% determine the top-level directory to use
%-------------------------------------------------------------------------------

[url topdir] = UFlocation ;

%-------------------------------------------------------------------------------
% export the matrices
%-------------------------------------------------------------------------------

for id = list

    % get the MATLAB version
    Problem = UFget (id, index) ;
    disp (Problem) ;

    % create the MM and RB versions
    UFwrite (Problem, [topdir 'MM'], 'MM', 'tar') ;
    UFwrite (Problem, [topdir 'RB'], 'RB', 'tar') ;

    % check the new MM and RB versions
    if (check)
	for format = { 'MM' , 'RB' }
	    try
		if (isempty (tmp))
		    P2 = UFread ([topdir format{1} filesep Problem.name]) ;
		else
		    P2 = UFread ([topdir format{1} filesep Problem.name], tmp) ;
		end
	    catch
		% The Problem may be too large for two copies to be in the
		% MATLAB workspace at the same time.  This is not an error,
		% but it means that the Problem cannot be checked.
		P2 = [ ] ;
		fprintf ('Unable to read %s/%s\n', format {1}, Problem.name) ;
		fprintf ('%s\n', lasterr) ;
	    end
	    if (~isempty (P2) && ~isequal (Problem, P2))
                Problem
                P2
		error ('%s version mismatch: %s\n', format {1}, Problem.name) ;
	    end
	    clear P2
	end
    end
end


text:
function [bn,burstiness,mburstf,mibsn,mbl,mibfr,mibif,burstf,ibsn,burstlen,ibfr,ibif,bp,normbnbdiff,burst,burstspikes,nonburstspikes] = burst_cls_kmeans(segspikes,sro)

% K-means clustering of ISIs
tivs = cell(1,size(segspikes,2));
ivnum = zeros(1,size(segspikes,2));
segfirstiv = zeros(1,size(segspikes,2));
seglastiv = zeros(1,size(segspikes,2));
for b1 = 1 : size(segspikes,2)
    tivs{b1} = diff(segspikes{b1});
    ivnum(b1) = length(tivs{b1});
    if b1 == 1
        seglastiv(b1) = ivnum(b1);
    else
        seglastiv(b1) = sum(ivnum(1 : b1 - 1)) + ivnum(b1);
    end
    if b1 == 1
        segfirstiv(b1) = 1;
    else
        segfirstiv(b1) = seglastiv(b1 - 1) + 1;
    end
end
ivs = cat(2,tivs{:});

%k-means starts
if length(ivs) > 49
    replino = 5;
    maxclust = 5;
    qual = zeros(1,maxclust);
    clustid = zeros(length(ivs),maxclust);
    for kmc1 = 2 : maxclust
        clustid(:,kmc1) = kmeans(ivs,kmc1,'emptyaction','drop','replicates',replino);
        qual(kmc1) = mean(silhouette(ivs',clustid(:,kmc1)));
    end
    ici = find(qual == max(qual));
    idealclust = clustid(:,ici);
    isiclusters = cell(1, ici);
    isiclustmean = zeros(1,ici);
    for kmc2 = 1 : ici
        isiclusters{kmc2} = ivs(idealclust == kmc2);
        isiclustmean(kmc2) = mean(isiclusters{kmc2});
    end
    ibclust = isiclusters{isiclustmean == min(isiclustmean)};
%k-means ends

%bf
    if isempty(ibclust) == 0
        firstspikes = cell(1,size(segspikes,2));
        lastspikes = cell(1,size(segspikes,2));
        segburstivs = cell(1,size(segspikes,2));
        for li1 = 1 : size(segspikes,2)
            segburstivs{li1} = find(tivs{li1} <= max(ibclust));
            sbiv1 = diff(segburstivs{li1});
            sbiv2 = find(sbiv1 > 1);
            firstspikes{li1} = segspikes{li1}(segburstivs{li1}(sbiv2 + 1));
            lastspikes{li1} = segspikes{li1}(segburstivs{li1}(sbiv2) + 1);
            firstspikes{li1} = [segspikes{li1}(1) firstspikes{li1}];
            lastspikes{li1} = [lastspikes{li1} segspikes{li1}(end)];
        end
        burst = [cat(2,firstspikes{:});cat(2,lastspikes{:})];
        burstlen = (burst(2,:) - burst(1,:)) / (sro / 1000); %length of bursts in ms
        burst(:,burstlen > 400) = [];
        burstlen(burstlen > 400) = [];
%bf

% maxthetacycle = 400;
% nonbursts = burst(:,bli >= maxthetacycle);
% burst(:,bli >= maxthetacycle) = [];

%Selection of burst and non-burst spikes
        spktm = cat(2,segspikes{:});
        burstspikes = cell(1,size(burst,2));
        ibsn = zeros(1,size(burst,2));
        for bj = 1 : size(burst,2)
            burstspikes{bj} = spktm(spktm >= burst(1,bj) & spktm <= burst(2,bj));
            ibsn(bj) = length(burstspikes{bj});
        end
        burstspkvect = cat(2,burstspikes{:});
        nbs = zeros(1,length(spktm));
        for nbi = 1 : length(spktm)
            for nbj = 1 : length(burstspkvect)
                if spktm(nbi) == burstspkvect(nbj)
                nbs(nbi) = 1;
                end
            end
        end
        nonburstspikes = spktm(nbs == 0);
%selection of burst & non-burst intervals
        nbss = zeros(1,length(ivs));
        for ni = 1 : length(ivs)
            for nj = 1 : length(ibclust)
                if ivs(ni) == ibclust(nj)
                    nbss(ni) = 1;
                end
            end
        end
        nbivs = ivs(nbss ~= 1);

%Comparison of burst and non-burst intervals
        if length(nbivs) >= 5
            bp = ranksum(ibclust,nbivs,0.01,'tail','both'); %ranksum test for burst - non-burst intervals
        else
            bp = -1;
        end
        normbnbdiff = (median(nbivs) - median(ibclust)) / (median(nbivs) + median(ibclust)); %normalized diff of burst - non-burst intervals

%burst parameter definitions
        mbl = median(burstlen); %median burst length
        ibif = 1./ibclust * sro; %intraburst inst.frequency
        mibif = median(1./ibclust*sro); %median intraburst inst.frequency
        ibfr = (ibsn ./ (burstlen/1000)); %intraburst firing rate per burst
        mibfr = median(ibfr); %median intraburst firing rate
        burstiness = (size(burstspkvect,2) / size(cat(2,segspikes{:},2),2)) * 100; %proportion of burst spikes
        covibsi = std(ibclust) / mean(ibclust); %COV of intraburst intervals
        first_isis = diff(cat(2,firstspikes{:})); %inter-firstspike intervals
        mfisis = median(first_isis) / (sro / 1000); % median interburst interval (ms)
        burstf = (1./first_isis)*sro; %inter-firstspike IF
        mburstf = median(burstf); %burst occurrence frequency
        bn = size(burst,2); %number of bursts
        mibsn = median(ibsn); %median intraburst spike number
    else
        bn = NaN;
        burstiness = NaN;
        mburstf = NaN;
        mibsn = NaN;
        mbl = NaN;
        mibfr = NaN;
        mibif = NaN;
        burstf = NaN;
        ibsn = NaN;
        burstlen = NaN;
        ibfr = NaN;
        ibif = NaN;
        bp = NaN;
        normbnbdiff = NaN;
        burst = NaN;
        burstspikes = NaN;
        nonburstspikes = NaN;
    end
else
    bn = NaN;
    burstiness = NaN;
    mburstf = NaN;
    mibsn = NaN;
    mbl = NaN;
    mibfr = NaN;
    mibif = NaN;
    burstf = NaN;
    ibsn = NaN;
    burstlen = NaN;
    ibfr = NaN;
    ibif = NaN;
    bp = NaN;
    normbnbdiff = NaN;
    burst = NaN;
    burstspikes = NaN;
    nonburstspikes = NaN;
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
filename='BridgeCool_Quadrilateral_Bilinear_Structured';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance'};
weights = 1;
constraint = {'volume'};
optimizer = 'MMA'; incrementFactor = 1;
filterType = 'P1';
constraint_case = 'INEQUALITY';
line_search_initiator = 'STANDARD';

nsteps = 15;
Vfrac_final = 0.2;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%% downloads

% data from: http://download.alleninstitute.org/informatics-archive/current-release/mouse_ccf/

% nrrd reader from: https://uk.mathworks.com/matlabcentral/fileexchange/34653-nrrd-format-file-reader

% structure tree from:
% http://api.brain-map.org/api/v2/data/query.csv?criteria=model::Structure,rma::criteria,[ontology_id$eq1],rma::options[order$eq%27structures.graph_order%27][num_rows$eqall]

%% sanitize structure tree

sanitizeStructureTree('query.csv', 'structure_tree_safe_2017.csv');
st = loadStructureTree('structure_tree_safe_2017.csv');

%% convert annotation volume to index version
% This cell can take a while, both for loading and doing the conversion

% av = nrrdread('annotation_10.nrrd');
% st = loadStructureTree('structure_tree_safe.csv');

avI = av;
ind = st.index;
id = st.id;
idS = sparse(double(id), ones(size(id)), double(ind), double(max(id)),1);
avI(av==0) = 997;
tic; avS = idS(avI); toc
avS = full(avS);
avI = reshape(avS, size(av));
avI = uint16(avI+1);

% running this for the 2017 version, I have to also add this line:
avI = permute(avI, [2 1 3]);

% writeNPY(avI, 'annotation_volume_10um_by_index.npy')

%% make colormap from structure_tree based on index

q = st.color_hex_triplet;

q(cellfun(@numel,q)==5) = {'019399'}; % special case where leading zero was evidently dropped
c1 = cellfun(@(x)hex2dec(x(1:2)), q, 'uni', false);
c2 = cellfun(@(x)hex2dec(x(3:4)), q, 'uni', false);
c3 = cellfun(@(x)hex2dec(x(5:6)), q, 'uni', false);
cmap = horzcat(vertcat(c1{:}),vertcat(c2{:}),vertcat(c3{:}))./255;

% save allen_ccf_colormap.mat cmap

author: cortex-lab

repo: allenCCF

save_path: github-repos/MATLAB/cortex-lab-allenCCF


text:
function netOutputs = fromCompiledNet(net)
%FROMCOMPILEDNET Decompiles a Net back into Layer objects
%   OUTPUTS = Layer.fromCompiledNet(NET) decompiles a compiled network (of
%   class Net), into their original Layer objects (i.e., a set of
%   recursively nested Layer objects).
%
%   Returns a cell array of Layer objects, each corresponding to an output
%   of the network. These can be composed with other layers, or compiled
%   into a Net object for training/evaluation.

% Copyright (C) 2017 Joao F. Henriques.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

  % call a local function, to work around an occasional Matlab crash that
  % happens only for static methods
  netOutputs = fromCompiledNet_(net);
end

function netOutputs = fromCompiledNet_(net)
  % use local copies so that any changes won't reflect on the original Net
  forward = net.forward ;
  backward = net.backward ;

  % first undo ReLU short-circuiting, since it creates cycles in the graph
  [forward, backward] = undoShortCircuit(net, forward, backward);

  % main decompilation code.

  % create a cell array 'layers' that will contain all created Layers,
  % indexed by their output vars. e.g., if layer L outputs var V, then
  % layers{V} = L.
  assert(~isempty(net.vars)) ;
  layers = cell(size(net.vars)) ;

  % create Input layers
  inputNames = fieldnames(net.inputs) ;

  for k = 1:numel(inputNames)
    var = net.inputs.(inputNames{k}) ;  % the var index
    layers{var} = Input('name', inputNames{k}, 'gpu', net.isGpuVar(var)) ;
  end

  % create Param layers
  for k = 1:numel(net.params)
    p = net.params(k) ;

    layer = Param('name', p.name, 'value', net.vars{p.var}, ...
      'gpu', net.isGpuVar(p.var), 'learningRate', p.learningRate, ...
      'weightDecay', p.weightDecay) ;

    layer.trainMethod = layer.trainMethods{p.trainMethod} ;  % index to string
    layers{p.var} = layer ;
  end

  % create the layers in forward-order. this order is important, because
  % when processing a layer all its inputs must have been processed too.
  for k = 1:numel(forward)
    layer = forward(k) ;
    args = layer.args ;

    % replace vars in the args list with the corresponding Layer objects
    for i = 1:numel(layer.inputVars)
      assert(~isempty(layers{layer.inputVars(i)}));  % assume already created
      args{layer.inputArgPos(i)} = layers{layer.inputVars(i)} ;
    end

    % now create a Layer with those arguments
    obj = Layer(layer.func, args{:}) ;
    obj.name = layer.name ;

    % retrieve the number of input derivatives from the backward struct
    obj.numInputDer = backward(end - k + 1).numInputDer ;

    %retrieve precious value (set during overloaded construction)
    obj.precious = layer.precious;

    % store in the 'layers' cell array so further references to the same
    % var will fetch the same layer
    obj.name = layer.name ;
    layers{layer.outputVar(1)} = obj ;

    % if the layer has multiple outputs, create a Selector for each of the
    % outputs after the first one
    for i = 2:numel(layer.outputVar)
      layers{layer.outputVar(i)} = Selector(obj, i) ;
    end
  end

  % the last var is the root (skipping the corresponding derivative var)
  rootLayer = layers{end - 1};

  if ~isequal(rootLayer.func, @root)
    % single output
    netOutputs = {rootLayer} ;
  else
    % root layer's arguments are the network's outputs
    netOutputs = rootLayer.inputs ;
  end

  % copy meta properties to the Layers
  for o = 1:numel(netOutputs)
    netOutputs{o}.meta = net.meta ;
  end
end

function [forward, backward] = undoShortCircuit(net, forward, backward)
  % undo ReLU short-circuiting, since it creates cycles in the graph.
  % collect unused vars left behind by short-circuiting, to reuse them now.
  unused = true(size(net.vars)) ;
  unused([forward.inputVars]) = false ;
  unused([forward.outputVar]) = false ;
  unused(2:2:end) = false ;  % ignore derivative vars
  unused = find(unused) ;  % convert to list of var indexes

  next = 1 ;  % next unused var to take
  for k = 1:numel(forward)
    if isequal(forward(k).func, @vl_nnrelu) && ...
      forward(k).inputVars(1) == forward(k).outputVar

      % a short-circuited ReLU (output var = input var). change its output
      % to a new unused var, and update any layers that depend on it.
      oldVar = forward(k).outputVar ;
      for j = k + 1 : numel(forward)
        in = forward(j).inputVars ;
        if any(in == oldVar)  % update this dependent layer
          forward(j).inputVars(in == oldVar) = unused(next) ;
        end
      end

      forward(k).outputVar = unused(next) ;  % update the ReLU

      next = next + 1 ;
      % shouldn't happen if short-circ. ReLUs leave unused vars behind
      assert(next <= numel(unused) + 1) ;
    end
  end
end

repo: MFT


text:
function neqn = p06_neqn ( )

%*****************************************************************************80
%
%% P06_NEQN returns the number of equations for problem P06.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer NEQN, the number of equations.
%
  neqn = 2;

  return
end


text:
Fid = 519;
Fnd = mesh.p(mesh.f(Fid,:),:);
idtmp = [1,2,3,1];
intpt0ID = mesh.eLoc(mesh.f_e(Fid,:));
intptID = -mesh.eLoc(mesh.f_e(Fid,intpt0ID<0));
intpt = mesh.eIntP(intptID,:);

plot3(Fnd(idtmp,1),Fnd(idtmp,2),Fnd(idtmp,3),'k','LineWidth',1)
hold on
plot3(intpt(:,1),intpt(:,2),intpt(:,3),'r','LineWidth',1)

testE1 = mesh.p(mesh.e(mesh.f_e(Fid,1),2),:) - intpt(1,:); % the first cutting edge corresponds to the first edge
mpt1 = 1/2*(mesh.p(mesh.e(mesh.f_e(Fid,1),2),:) + intpt(1,:));
testE2 = intpt(1,:) - mesh.p(mesh.e(mesh.f_e(Fid,1),1),:);
mpt2 = 1/2*(intpt(1,:) + mesh.p(mesh.e(mesh.f_e(Fid,1),1),:));
lE = norm(testE1) + norm(testE2);
ftId = -mesh.tLoc(mesh.f_t(Fid,:));
bas11 = squeeze(femI.bas1(ftId(1),:,:));
bas12 = squeeze(femI.bas2(ftId(1),:,:));
bas21 = squeeze(femI.bas1(ftId(2),:,:));
bas22 = squeeze(femI.bas2(ftId(2),:,:));

% use this plot to check the piece
% scatter3(mpt1(:,1),mpt1(:,2),mpt1(:,3),'r')

dof1 = zeros(6,1);
for i = 1:6

    dof1(i) = (dot((cross(bas11(1:3,i),mpt2') + bas11(4:6,i)),testE2)+ ...
        dot((cross(bas12(1:3,i),mpt1') + bas12(4:6,i)),testE1));

end

dof2 = zeros(6,1);
for i = 1:6

    dof2(i) = (dot((cross(bas21(1:3,i),mpt2') + bas21(4:6,i)),testE2)+ ...
        dot((cross(bas22(1:3,i),mpt1') + bas22(4:6,i)),testE1));

end
% can check which this edge is associated which dof, and which piece is
% related to bas1 or bas2
i1 = find(abs(dof1-1)<10^(-10)); i2 = find(abs(dof2-1)<10^(-10));
bas11 = bas11(:,i1); bas12 = bas12(:,i1);
bas21 = bas21(:,i2); bas22 = bas22(:,i2);

xm = min(Fnd(idtmp,1));
xM = max(Fnd(idtmp,1));
ym = min(Fnd(idtmp,2));
yM = max(Fnd(idtmp,2));
h = xM-xm;
[x,y] = meshgrid(xm:h/16:xM,ym:h/16:yM);
z = ones(size(x))*Fnd(1,3); % assume this face is parallel to xy plane
szh = size(x,1);
funx1 = zeros(size(x));
funy1 = zeros(size(x));
funz1 = zeros(size(x));

for i = 1:szh
    for j = 1:szh

        if x(i,j) < intpt(1,1)
            piece_ind = 1;
        elseif x(i,j) >= intpt(1,1)
            piece_ind = 2;
        end

        vec = ifebas(x(i,j), y(i,j),z(i,j),bas11,bas12,piece_ind);

        funx1(i,j) = vec(1);
        funy1(i,j) = vec(2);
        funz1(i,j) = vec(3);

        if y(i,j)-x(i,j)>Fnd(1,2)-Fnd(1,1)
            funx1(i,j) = 0; funy1(i,j) = 0; funz1(i,j) = 0;
        end

    end
end
quiver3(x,y,z+0.03,funx1,funy1,funz1,2,'b')

funx2 = zeros(size(x));
funy2 = zeros(size(x));
funz2 = zeros(size(x));

for i = 1:szh
    for j = 1:szh

        if x(i,j) < intpt(1,1)
            piece_ind = 1;
        elseif x(i,j) >= intpt(1,1)
            piece_ind = 2;
        end

        vec = ifebas(x(i,j), y(i,j),z(i,j),bas21,bas22,piece_ind);

        funx2(i,j) = vec(1);
        funy2(i,j) = vec(2);
        funz2(i,j) = vec(3);

        if y(i,j)-x(i,j)>Fnd(1,2)-Fnd(1,1)
            funx2(i,j) = 0; funy2(i,j) = 0; funz2(i,j) = 0;
        end

    end
end
quiver3(x,y,z-0.03,funx2,funy2,funz2,2,'color',[1,0.5,0])

box on
axis equal

Elem1Nd = mesh.p(mesh.t(mesh.f_t(Fid,1),:),:);
Elem2Nd = mesh.p(mesh.t(mesh.f_t(Fid,2),:),:);
idtmp = [1,2,3,1,4,3,1,4,2];
plot3(Elem1Nd(idtmp,1),Elem1Nd(idtmp,2),Elem1Nd(idtmp,3),'k','LineWidth',1)
hold on
plot3(Elem2Nd(idtmp,1),Elem2Nd(idtmp,2),Elem2Nd(idtmp,3),'k','LineWidth',1)
Eintpt1Id = mesh.eLoc(mesh.t_e(mesh.f_t(Fid,1),:));
Eintpt1Id = -mesh.eLoc(mesh.t_e(mesh.f_t(Fid,1),Eintpt1Id<0));
Eintpt1 = mesh.eIntP(Eintpt1Id,:);
Eintpt2Id = mesh.eLoc(mesh.t_e(mesh.f_t(Fid,2),:));
Eintpt2Id = -mesh.eLoc(mesh.t_e(mesh.f_t(Fid,2),Eintpt2Id<0));
Eintpt2 = mesh.eIntP(Eintpt2Id,:);
patch(Eintpt1(:,1),Eintpt1(:,2),Eintpt1(:,3),'r','FaceAlpha',0.7)
idtmp = [1,2,4,3];
patch(Eintpt2(idtmp,1),Eintpt2(idtmp,2),Eintpt2(idtmp,3),'r','FaceAlpha',0.7)
box on
axis off
axis equal

function u = ifebas(x,y,z,bas1,bas2,piecie)
fnorm = [0;0;1]; % assume this face is parallel to xy plane

if piecie == 1
    bas = bas1;
elseif piecie == 2
    bas = bas2;
end

u = cross(bas(1:3),[x;y;z]) + bas(4:6);
u = cross(u,fnorm);

end

repo: ifem


text:
function denoised_im = denoise_foe(im, basis, W, a, sigma)
%DENOISE_FOE Denoise image using Field of Experts prior model
%
%   denoised_im = denoise_foe(im, basis, W, a, sigma)
%
% Denoises the image im, using a Gaussian noise model with variance sigma,
% and a Field of Experts image prior specified by basis, W, and a.
%
%
% (C) Laurens van der Maaten, 2009
% Delft University of Technology

    % Initialize some variables
    eta = 10;
    max_iter = 275; %275 je bilo originalno
    im = double(im);
    dI = zeros(size(im));

    % Convert to YCbCr color space if input is RGB
    if size(im, 3) == 3
        im = 255 .* rgb2ycbcr(im ./ 255);
    end
    denoised_im = im;

    minEnergy = realmax;
    minEnergyImg = denoised_im;
    % Perform gradient iterations
    for iter=1:max_iter

        % Compute gradient of posterior log-likelihood
        for c=1:size(denoised_im, 3)
            dI(:,:,c) = reshape(foe_energy_grad_x(denoised_im(:,:,c), basis, W, a), [size(im, 1) size(im, 2)]);
        end
        dI = dI + (1 / sigma ^ 2) * (im - denoised_im);

        % Perform gradient update
        denoised_im = denoised_im - eta * dI;

        % Print progress
        if ~rem(iter, 10) || iter == max_iter

            % Compute energy of posterior
            E = foe_energy(denoised_im, basis, W, a, size(denoised_im)) + (1 ./ (sigma ^ 2)) * sum(((im(:) - denoised_im(:)) ./ 255) .^ 2);

            if E < minEnergy
                minEnergy = E;
                minEnergyImg = denoised_im;
                disp(['Iteration ' num2str(iter) ': energy is ' num2str(E)]);
            else
                disp('Finished because of energy increase.');
                denoised_im = minEnergyImg;
                break;
            end

            % Show intermediate result
            if size(im, 3) == 3
                tmp1 = uint8(round(255 .* ycbcr2rgb(im ./ 255)));
                tmp2 = uint8(round(255 .* ycbcr2rgb(denoised_im ./ 255)));
            else
                tmp1 = uint8(round(im));
                tmp2 = uint8(round(denoised_im));
            end
        end
    end

    subplot(1, 2, 1); imshow(tmp1); colormap(gray); title('Noisy image');
    subplot(1, 2, 2); imshow(tmp2); colormap(gray); title(['Denoised image (iteration ' num2str(iter) ')']);
    drawnow;

    % Convert back to RGB if necessary
    if size(denoised_im, 3) == 3
        denoised_im = 255 .* ycbcr2rgb(denoised_im ./ 255);
    end
    denoised_im = uint8(round(denoised_im));

author: lbasek

repo: image-denoising-benchmark


text:
classdef mme_bus_opf_acc < mp.mme_bus_opf_ac

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

%     properties
%         name = 'bus';
%     end

    methods
        function add_constraints(obj, mm, nm, dm, mpopt)
            dme = obj.data_model_element(dm);
            nme = obj.network_model_element(nm);

            %% voltage angle reference constraint
            ref = find(nme.node_types(nm, dm) == mp.NODE_TYPE.REF);
            varef = dme.va_start(ref);
            fcn_vref = @(xx)va_fcn(nme, xx, ref, varef);
            hess_vref = @(xx, lam)va_hess(nme, xx, lam, ref);
            mm.add_nln_constraint('Vref', length(ref), 1, fcn_vref, hess_vref, {'Vr', 'Vi'});

            %% fixed voltage magnitudes
            veq = find(dme.vm_lb == dme.vm_ub);
            nveq = length(veq);
            if nveq
                fcn_vm2eq = @(xx)vm2_fcn(nme, xx, veq, dme.vm_ub(veq).^2);
                hess_vm2eq = @(xx, lam)vm2_hess(nme, xx, lam, veq);
                mm.add_nln_constraint('Veq', nveq, 1, fcn_vm2eq, hess_vm2eq, {'Vr', 'Vi'});
            end
            mm.userdata.veq = veq;

            %% voltage magnitude limits
            viq = find(dme.vm_lb ~= dme.vm_ub);
            nviq = length(viq);
            if nviq
                fcn_vlim = @(xx)vm2_fcn(nme, xx, viq, ...
                        {dme.vm_lb(viq).^2, dme.vm_ub(viq).^2} );
                hess_vlim = @(xx, lam)vm2_hess(nme, xx, lam, viq);
                mm.add_nln_constraint({'Vmin', 'Vmax'}, [nviq;nviq], 0, fcn_vlim, hess_vlim, {'Vr', 'Vi'});
            end
            mm.userdata.viq = viq;
        end

        function x0 = interior_x0(obj, mm, nm, dm, x0)
            vv = mm.get_idx();
            varef1 = mm.interior_va(nm, dm);
            vm = obj.interior_vm(mm, nm, dm);
            v_ = vm * exp(1j*varef1);
            x0(vv.i1.Vr:vv.iN.Vr) = real(v_);
            x0(vv.i1.Vi:vv.iN.Vi) = imag(v_);
        end

        function obj = data_model_update(obj, mm, nm, dm, mpopt)
            %% complex bus voltages
            nn = nm.get_idx('node');
            V = nm.soln.v(nn.i1.bus:nn.iN.bus);

            %% shadow prices on voltage magnitudes
            [nne, nni] = mm.get_idx('nle', 'nli');
            lambda = mm.soln.lambda;
            mu_vm_lb = zeros(nn.N.bus, 1);  %% init to all 0
            mu_vm_ub = mu_vm_lb;            %% init to all 0
            if mm.userdata.veq
                lam = lambda.eqnonlin(nne.i1.Veq:nne.iN.Veq);
                lam_p = zeros(size(lam));
                lam_n = zeros(size(lam));
                lam_p(lam > 0) =  lam(lam > 0);
                lam_n(lam < 0) = -lam(lam < 0);
                mu_vm_lb(mm.userdata.veq) = lam_n;
                mu_vm_ub(mm.userdata.veq) = lam_p;
            end
            mu_vm_lb(mm.userdata.viq) = lambda.ineqnonlin(nni.i1.Vmin:nni.iN.Vmin);
            mu_vm_ub(mm.userdata.viq) = lambda.ineqnonlin(nni.i1.Vmax:nni.iN.Vmax);

            vm = abs(V);
            mu_vm_lb = mu_vm_lb .* vm * 2;
            mu_vm_ub = mu_vm_ub .* vm * 2;

            %% shadow prices on node power balance
            [lam_p, lam_q] = mm.node_power_balance_prices(nm);
            lam_p = lam_p(nn.i1.bus:nn.iN.bus);     %% for bus nodes only
            lam_q = lam_q(nn.i1.bus:nn.iN.bus);     %% for bus nodes only

            %% update in the data model
            dme = obj.data_model_element(dm);
            dme.tab.va(dme.on) = angle(V) * 180/pi;
            dme.tab.vm(dme.on) = abs(V);
            dme.tab.lam_p(dme.on) = lam_p / dm.base_mva;
            dme.tab.lam_q(dme.on) = lam_q / dm.base_mva;
            dme.tab.mu_vm_lb(dme.on) = mu_vm_lb;
            dme.tab.mu_vm_ub(dme.on) = mu_vm_ub;
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower