text:
function [Xco Tau MT] = HC_CrankNicolson_h (Tin,TC,DE,HC,SE,DI,TS,TT,TBC,BC)
%
% DESCRIPTION:
%
% Crank-Nicolson numerical method for one dimensional unsteady state
% heat transfer by conduction for homogenous material
%
% Basic PDE equation :
%
%     dT         d^2 T
%   ------ = TD ------             for 0 < x < DI  and  tau > 0
%    dtau        dx^2
%
% Principle of method :
%       ^
%  time |
%       |      m-1     m     m+1
%       |       |      |      |
%       |    ---o------x------o-----   p+1
%       |  dtau |      |      |
%       |    ---o------o------o-----   p
%       |  dtau |      |      |
%       |    ---+------+------+-----   p-1
%       |       |  dx  |  dx  |
%       |
%       +-------------------------------->
%                                        x
% Space description:
%                  DI
%   |<------------------------------>|
%   +   +   +   +   +  . . .     +   +
%   1   2   3   4   5           n-1  n
%   |<->|
%     SE
%
% Difference equation :
%
%   T(m,p+1) - T(m,p)     TD
%  ------------------- = -------( T(m-1,p+1) - 2.T(m,p+1) + T(m+1,p+1) )
%       dTau             2.dx.dx
%
%                         TD
%                      + -------( T(m-1,p) - 2.T(m,p) + T(m+1,p) )
%                        2.dx.dx
%
%  System of equations:
%
% (1+M).T(2,p+1) - M/2.T(3,p+1) =
%                 =  M/2.T(1,p) + (1-M).T(2,p) + M/2.T(3,p) +  M/2.T(1,p+1)
%  . . .
%  for m <2,n-1>
% - M/2.T(m-1,p+1) + (1+M).T(m,p+1)- M/2.T(m+1,p+1) =
%                             =  M/2.T(m-1,p) + (1-M).T(m,p) + M/2.T(m+1,p)
%  . . .
%  - M/2.T(n-2,p+1) + (1+M).T(n-1,p+1) =
%              =  M/2.T(n-2,p) + (1-M).T(n-1,p) + M/2.T(n,p) + M/2.T(n,p+1)
%
%  Matrix notation:
%
%   A.T = b
%
%  The elements of A are given by
%
%   A(n-2,n-2) = 0 except for:
%
%   A(1,1)     = 1 + M;
%   A(1,2)     = - M/2;
%
%   A(m,m-1)   = - M/2;
%   A(m,m)     = 1 + M;
%   A(m,m+1)   = - M/2;   for m <2,n-3>
%
%   A(n-2,n-3) = - M/2;
%   A(n-2,n-2) = 1 + M;
%
%  and the elements of b are given by
%
%   b(1)   = M/2*Tin(1)   + (1-M)*Tin(2)   + M/2*Tin(3) + M/2*T(1);
%
%   b(m)   = M/2*Tin(m)   + (1-M)*Tin(m+1) + M/2*Tin(m+2);    for m <2,n-3>
%
%   b(n-2) = M/2*Tin(n-2) + (1-M)*Tin(n-1) + M/2*Tin(n) + M/2*T(n);
%
% INPUTS:
%
% Tin  - initialization temperatures                        [ K           ]
% TC   - thermal conductivity                               [ W.m^-1.K^-1 ]
% DE   - density                                            [ kg.m^-3     ]
% HC   - specific heat capacity                             [ kg.m^-3     ]
% SE   - size of element                                    [ m           ]
% DI   - distance                                           [ m           ]
% TS   - time step                                          [ s           ]
% TT   - total time of simulation                           [ s           ]
% TBC  - type of boundary condition
%        (1 = first-type, 2 = second-type)                  [ -           ]
% BC   - boundary condition
%         - first-type : BC(1) = T(1),  BC(2) = T(n)        [ K           ]
%         - second-type: BC(1) = iQ(1), BC(2) = iQ(n)       [ W.m^-2      ]
%
% OUTPUTS:
%
% Xco  - vector of X coordinate                             [ m           ]
% Tau  - vector of time                                     [ s           ]
% MT   - matrix of temperatures                             [ K           ]
%
% AUXILIARY VARIABLE:
%
% r    - time  step  number                                 [ -           ]
% n    - space point number                                 [ -           ]
% TD   - thermal diffusivity                                [ m^2.m^-1    ]
% M    - modul (Fourier number)                             [ -           ]
% A    - matrix (n-2,n-2)                                   [ -           ]
% b    - vector (n-2)                                       [ -           ]
%
% Copyright (C) 2013, Technical University of Kosice
% Author  : Zecova Monika, Terpak Jan
% Revision: 27.08.2013
%
  r = round(TT/TS) + 1;
  n = round(DI/SE) + 1;

  Xco = 0.0:SE:DI;
  Tau = 0.0:TS:TT;
  MT  = zeros(r,n);
  MT(1,:) = Tin;

  for j=2:r
    [T]     = CrankNicolson (Tin,TC,DE,HC,SE,TS,TBC,BC);
    Tin     = T;
    MT(j,:) = Tin;
  end
end
function [T] = CrankNicolson (Tin,TC,DE,HC,SE,TS,TBC,BC)

  TD = TC/(DE*HC);
  M  = TD*TS/(SE*SE);

  T = Tin;

  if TBC == 1
     T(1)   = BC(1);
     T(end) = BC(2);
  end

  if TBC == 2
     T(1)   = T(2)     - BC(1)*SE/TC;
     T(end) = T(end-1) - BC(2)*SE/TC;
  end

  k = length(Tin)-2;
  A = zeros(k,k);
  b = zeros(k,1);

  for i=1:k
      b(i) = M/2*Tin(i) + (1-M)*Tin(i+1)   + M/2*Tin(i+2) ;
  end

  A(1,1) = (1+M);
  A(1,2) = - M/2;
  b(1)   = b(1) + M/2*T(1);

  for i=2:k-1
      A(i,i-1) = - M/2;
      A(i,i)   = 1 + M;
      A(i,i+1) = - M/2;
  end

  A(k,k-1) = - M/2;
  A(k,k)   = 1 + M;
  b(k)     = b(k) + M/2*T(end);

  x = A\b;

  for i=2:length(Tin)-1
      T(i) = x(i-1);
  end

end


text:
%%Generate a simualted trajectory
%rstat: empty --> randomize (also return the seed in case it may be needed for debugging)
%       else --> use rstat as the seed (must be either a postive int or an rstat returned by rsg)
%dir_name : Directory to write outputs (must include trailing slash)
%ini_pva :[pos_in_nav(lLh), vel_in_body, body_nav_euler];
%out_typ :[[simulation_freq_multiplier imu_output_freq];[1 gps_freq]; [2 odo_freq]]
%sim_mode :
%   0:nav frame simulation,
%   1:simulation in a planar frame with constant g (position output in meters) (ini_pos must also be in meters)
%mot_def : Definition of motions. Each row defines a segment.
%   mot_def(:,1)=motion type (see the code for the options)
%   mot_def(:,2:5)=parameters of the motion (see the code)
%   mot_def(:,6)=maximum time for the given segment
%sen_conf: configuration for additional individual sensors (senconf=[la_in_master;euler_in_master])
%imu_conf: configuration for additional imus (imuconf=[la_in_master;euler_in_master])

%Output Files:
%mimu.bin: imu_data=[out_ind, acc, rotation_rate] (time=out_ind/imu_freq)
%mnav.bin: true_nav_data=[out_ind, pos_n or pos_dn, vel_b;att_n] (sim_mode determines position)
function rstat=PathGen(dir_name, ini_pva, mot_def, out_typ, sim_mode, sen_conf, imu_conf, vib_param, rstat)

if (rstat(1)~=0)
    randn('state',rstat);
else
    rstat=randn('state');
end

DEBUG=0;    %runs an additional conventional ins to see how compatible the trajectory generator and ins outputs are. (usually required to see the effect numerical approximations under vibration)

%% Frequencies
outfreq=out_typ(1,2);    %imu output frequency
intmul=out_typ(1,1);      %simulation freq multiplier. simulation freq=outfreq*intmul
dt=1/outfreq/intmul;

%%Additional Sensor configurations
nsen=0;     %number of additional sensor units
if (~isempty(sen_conf))    %there exist additional sensors
    nsen=size(sen_conf,1);
    sen_la=zeros(3,nsen);
    sen_or=zeros(3,3,nsen); %dcm from master to sensor
    sen_typ=sen_conf(:,7);
    for in=1:nsen
        sen_la(:,in)=sen_conf(in,1:3); %lever arm in the master body frame
        sen_or(:,:,in)=euler2dcm_v000(sen_conf(in,4:6)); %transformation matrix from sensor to master (euler angles defined in master);
    end
end

%%Additional imu configurations;
nimu=0;     %number of additional imu units
if (~isempty(imu_conf))    %there exist additional sensors
    nimu=size(imu_conf,1);
    imu_la=zeros(3,nimu);
    imu_or=zeros(3,3,nimu); %dcm from master to sensor
    imu_vib=zeros(3,3,nimu);
    for in=1:nimu
        imu_la(:,in)=imu_conf(in,1:3); %lever arm in the master body frame
        imu_or(:,:,in)=euler2dcm_v000(imu_conf(in,4:6)); %transformation matrix from sensor to master (euler angles defined in master);
        imu_vib(:,:,in)=eye(3);     %default value=no vibration
    end
end

%%vibration parameters for additional imus
if (nimu>0)
    vib_state=zeros(3,nimu);
    vib_model=zeros(2,nimu);
    if (~isempty(vib_param))
        for in=1:nimu
            vib_model(:,in)=vib_param(in,1:2)';
            vib_state(3,in)=randn(1)*vib_param(in,3);     %initial vibration angles
            imu_vib(:,:,in)=euler2dcm_v000(vib_state(:,in));
        end
    end
end

%% Path Gen Command Filter and feedback parameters (somehow arbitrary)
INP_FILT_A=eye(3)*0.9;
INP_FILT_B=eye(3)-INP_FILT_A;
att_n_dot=zeros(3,1);
vel_b_dot=zeros(3,1);

%% Open output files
fmimu=fopen([dir_name 'mimu.bin'],'wb');  %sensor outputs
fmnav=fopen([dir_name 'mnav.bin'],'wb');  %master navigation solution

if (nimu>0) %for each additional imus, create additional navigation and imu files
    fsnav=zeros(nimu,1);    %navigation solution files (for each full imu)
    fsimu=zeros(nimu,1);    %imu files
    for in=1:nimu
        fsnav(in)=fopen([dir_name 'snav' int2str(in) '.bin'],'wb');
        fsimu(in)=fopen([dir_name 'simu' int2str(in) '.bin'],'wb');
    end
end

nout=size(out_typ,1);
for in=2:nout    %measurement files
    switch out_typ(in,1)
        case {1}
            fgps=fopen([dir_name 'gps.bin'],'wb');
        case {2}
            fodo=fopen([dir_name 'odo.bin'],'wb');
    end
end

%% convert times into indeces
for in=2:nout;
    out_typ(in,2)=intmul*round(outfreq/out_typ(in,2));
end
for in=1:size(mot_def,1)
    mot_def(in,6)=round(mot_def(in,6)*outfreq*intmul);
end

%%%%%%%%%%% Start Computations %%%%%%%%%%%%%%%
sim_cnt=0;
fflag=1;
acc_tot=zeros(3,1);
gyro_tot=zeros(3,1);
odo_dist=0;

%master navigator states
att_n=ini_pva(:,3);
Cbn=euler2dcm_v000(att_n);
pos_n=ini_pva(:,1);
pos_dn=zeros(3,1);  %total position change
vel_b=ini_pva(:,2);
vel_n=Cbn*vel_b;
[Rn, Re, g, sL, cL, WIE]=geoparam_v000(pos_n);

%Write initial values to the file
fwrite(fmimu,[0;zeros(3,1);zeros(3,1);zeros(nsen,1)],'double');
if (sim_mode==0)
    vr_a=[0;pos_n+pos_dn;vel_b;att_n];
else
    vr_a=[0;pos_dn;vel_b;att_n];
end
fwrite(fmnav, vr_a,'double');

if (nimu>0)
    if (sim_mode==0)
        [spos_n, svel_b, sCsn, simu]=la_imu_v000(zeros(3,1), zeros(3,1), zeros(3,1), pos_n, vel_b, Cbn, imu_la, imu_or, 0, zeros(3,1), imu_vib);
    elseif (sim_mode==1)
        [spos_n, svel_b, sCsn, simu]=la_imu_v000(zeros(3,1), zeros(3,1), zeros(3,1), pos_dn, vel_b, Cbn, imu_la, imu_or, 1, zeros(3,1), imu_vib);
    end
    for in1=1:nimu
        satt_n=dcm2euler_v000(sCsn(:,:,nimu));
        fwrite(fsnav(in1), [0;spos_n(:,in1);svel_b(:,in1); satt_n], 'double');
        fwrite(fsimu(in1), [0;simu], 'double');
    end
end

%%variables for debugging
if (DEBUG)
    %strapdown states (mostly for debugging purposes)
    strap_Cbn=Cbn;
    strap_vel_b=vel_b;
    if (sim_mode==0)
        strap_pos_n=pos_n;
    elseif (sim_mode==1)
        strap_pos_n=pos_dn;
    end

    if (nimu>0)
        strap_spos_n=spos_n;
%         for in=1:nimu       %!!!THIS DOES DOT SEEM OK.!!!!
%             strap_svel_b=svel_b(:,in);
%         end
        strap_svel_b=svel_b;
        strap_sCsn=sCsn;
    end

    sr_a=sum(mot_def(:,6)); %maximum number of outputs
    deb_vel=zeros(3,sr_a);
    deb_pos=zeros(3,sr_a);
    deb_att=zeros(3,sr_a);
    deb_ctr=0;
end

%Start generation of trajectory
for in=1:size(mot_def,1)
    mmod=mot_def(in,1);

    if (mmod==1)
        %Determine the inputs directly
        att_n_dot_com=[mot_def(in,2);mot_def(in,3);mot_def(in,4)];
        vel_b_dot_com=[mot_def(in,5);0;0];
    elseif mmod==2    %new abs. att and vel values to reach
        %determine the commands
        att_n_com=[mot_def(in,2);mot_def(in,3);mot_def(in,4)];
        vel_b_com=[mot_def(in,5);0;0];
    elseif mmod==3    %rel att and vel changes
        att_n_com=[mot_def(in,2)+att_n(1);mot_def(in,3)+att_n(2);mot_def(in,4)+att_n(3)];
        vel_b_com=vel_b+[mot_def(in,5);0;0];
    elseif mmod==4    %abs attitude, rel vel
        att_n_com=[mot_def(in,2);mot_def(in,3);mot_def(in,4)];
        vel_b_com=vel_b+[mot_def(in,5);0;0];
    elseif mmod==5    %rel att, abs vel
        att_n_com=[mot_def(in,2)+att_n(1);mot_def(in,3)+att_n(2);mot_def(in,4)+att_n(3)];
        vel_b_com=[mot_def(in,5);0;0];
    end

    att_n_com_filt=att_n;
    vel_b_com_filt=vel_b;

    seg_cnt_lim=sim_cnt+mot_def(in,6);
    seg_next=0;

    while ((sim_cnt<seg_cnt_lim) && ~seg_next)
        %determine the inputs
        if (mmod==1)
            %filter the input
            att_n_dot=INP_FILT_A*att_n_dot+INP_FILT_B*att_n_dot_com;
            vel_b_dot=INP_FILT_A*vel_b_dot+INP_FILT_B*vel_b_dot_com;
        else
            %filter the command
            att_n_com_filt=INP_FILT_A*att_n_com_filt+INP_FILT_B*att_n_com;
            vel_b_com_filt=INP_FILT_A*vel_b_com_filt+INP_FILT_B*vel_b_com;

            %close the loop
            att_n_dot=att_n_com_filt-att_n;
            vel_b_dot=vel_b_com_filt-vel_b;

            if (norm(att_n_dot)<1e-4 && norm(vel_b_dot)<1e-4)
                seg_next=1;
            end
        end

        %Compute the master imu solution
        [acc, gyro, vel_n_dot, pos_n_dot, Cbn_dot]=comp_master(pos_n+pos_dn, vel_b, att_n, Cbn, vel_b_dot, att_n_dot, sim_mode, g);

        %update master results
        acc_tot=acc_tot+acc;
        gyro_tot=gyro_tot+gyro;
        pos_dn=pos_dn+pos_n_dot*dt;
        att_n=att_n+att_n_dot*dt;
        %Cbn=euler2dcm_v000(att_n);
        Cbn=Cbn+Cbn_dot*dt;
        vel_b=vel_b+vel_b_dot*dt;
        vel_n=Cbn*vel_b;
%         vel_n=vel_n+vel_n_dot*dt;
%         vel_b=Cbn'*vel_n;

        %%odometer resuts
        odo_vel=vel_b;
        odo_dist=odo_dist+norm(vel_b)*dt;

        %update simulation counter
        sim_cnt=sim_cnt+1;

        %write the results
        if (mod(sim_cnt,intmul)==0)
            %%average gyro and acc for the output interval
            gyro_avg=gyro_tot/intmul;
            acc_avg=acc_tot/intmul;

            if (fflag)  %only once: set gyro derivative to zero, write initial values to the files
                fflag=0;
                gyro_pre=gyro_avg;
                acc_pre=acc_avg;
            end

            %%Compute the sensor outputs
            gyro_der=(gyro_avg-gyro_pre)*outfreq;
            if (nsen>0)
                sensor_out=la_sen_v000(acc_avg, gyro_avg, gyro_der, sen_la, sen_or, sen_typ);
            else
                sensor_out=[];
            end
            fwrite(fmimu,[(sim_cnt/intmul);acc_avg;gyro_avg;sensor_out],'double');
            if (sim_mode==0)
                vr_a=[(sim_cnt/intmul);pos_n+pos_dn;vel_b;att_n];
            else
                vr_a=[(sim_cnt/intmul);pos_dn;vel_b;att_n];
            end
            fwrite(fmnav, vr_a,'double');

            %Compute the slave imu outputs
            if (nimu>0)
                %update the vibration
                [gyro_vib, vib_state, imu_vib_new]=update_vib(vib_state, vib_model, imu_vib, imu_or, 1/outfreq);

                %compute the slave results
                if (sim_mode==0)
                    [spos_n, svel_b, sCsn, simu]=la_imu_v000(acc_avg, gyro_avg, gyro_der, pos_n+pos_dn, vel_b, Cbn, imu_la, imu_or, 0, gyro_vib, imu_vib);
                elseif (sim_mode==1)
                    [spos_n, svel_b, sCsn, simu]=la_imu_v000(acc_avg, gyro_avg, gyro_der, pos_dn, vel_b, Cbn, imu_la, imu_or, 1, gyro_vib, imu_vib);
                end
                for in1=1:nimu
                    satt_n=dcm2euler_v000(sCsn(:,:,nimu));
                    fwrite(fsnav(in1), [(sim_cnt/intmul);spos_n(:,in1);svel_b(:,in1);satt_n(:,in1)], 'double');
                    fwrite(fsimu(in1), [(sim_cnt/intmul);simu(:,in1)], 'double');
                end
                imu_vib=imu_vib_new;
            end

            %prepare for the next cycle
            gyro_pre=gyro_avg;
            gyro_tot=zeros(3,1);
            acc_tot=zeros(3,1);

            %%debug
            if (DEBUG)
                %%run the sptradown
                if (sim_mode==0)
                    [strap_Cbn, strap_vel_b, strap_pos_n]=strapdown_bd_dcm_v000(strap_Cbn, strap_vel_b, strap_pos_n, acc_avg, gyro_avg, 1/outfreq);
                    for in1=1:nimu
                        [strap_sCsn(:,:,in1), strap_svel_b(:,in1), strap_spos_n(:,in1)]=strapdown_bd_dcm_v000(strap_sCsn(:,:,in1), strap_svel_b(:,in1), strap_spos_n(:,in1), simu(1:3,in1), simu(4:6,in1), 1/outfreq);
                    end
                elseif (sim_mode==1)
                    [strap_Cbn, strap_vel_b, strap_pos_n]=strapdown_pln_dcm_v000(strap_Cbn, strap_vel_b, strap_pos_n, acc_avg, gyro_avg, g, 1/outfreq, 1);
                    for in1=1:nimu
                        [strap_sCsn(:,:,in1), strap_svel_b(:,in1), strap_spos_n(:,in1)]=strapdown_pln_dcm_v000(strap_sCsn(:,:,in1), strap_svel_b(:,in1), strap_spos_n(:,in1), simu(1:3,in1), simu(4:6,in1),g, 1/outfreq,1);
                    end
                end

                deb_ctr=deb_ctr+1;
%                 deb_vel(:,deb_ctr)=strap_svel_b-svel_b;
%                 deb_pos(:,deb_ctr)=strap_spos_n-spos_n;
%                 deb_att(:,deb_ctr)=dcm2euler_v000(strap_sCsn(:,:,1)*sCsn');
                deb_vel(:,deb_ctr)=strap_vel_b-vel_b;
                deb_pos(:,deb_ctr)=strap_pos_n-pos_dn-pos_n;
                deb_att(:,deb_ctr)=dcm2euler_v000(strap_Cbn(:,:,1)*Cbn');
            end

        end

        %%measurement files
        if (nout>1)
            for sr_a=2:nout;
                switch out_typ(sr_a,1)
                    case {1}
                        if (mod(sim_cnt,out_typ(sr_a,2))==0)
                            if (sim_mode==0)
                                vr_a=[(sim_cnt/intmul);vel_n;pos_n+pos_dn];
                            elseif(sim_mode==1)
                                vr_a=[(sim_cnt/intmul);vel_n;pos_dn];
                            end
                            fwrite(fgps,vr_a,'double');
                        end
                    case {2}
                        if (mod(sim_cnt,out_typ(sr_a,2))==0)
                            fwrite(fodo, [(sim_cnt/intmul);odo_vel;odo_dist],'double');
                        end
                end
            end
        end
    end
end
if (DEBUG)
    deb_vel(:,deb_ctr)
    deb_pos(:,deb_ctr)
    deb_att(:,deb_ctr)
end
fclose('all');

%vel_typ==1 -> b, vel_typ==2 -> n
function [acc, gyro, vel_n_dot, pos_n_dot, Cbn_dot]=comp_master(pos_n, vel_b, att_n, Cbn, vel_b_dot, att_n_dot, styp, g)
%velocity in n
vel_n=Cbn*vel_b;

%%Calculate Spherical Earth Parameters and gravity
if (styp==0)
    [Rn, Re, g, sL, cL, WIE]=geoparam_v000(pos_n);
    Rn_eff=Rn+pos_n(3);
    Re_eff=Re+pos_n(3);
    gravity=[0;0;-g];
elseif (styp==1)
    gravity=[0;0;-g];
end

%compute w_nb_n (rotation rate of "b" wrt "n" defined in "n")
sh=sin(att_n(3));ch=cos(att_n(3));
w_nb_n=zeros(3,1);
w_nb_n(1)=-sh*att_n_dot(2)+Cbn(1,1)*att_n_dot(1);
w_nb_n(2)=ch*att_n_dot(2)+Cbn(2,1)*att_n_dot(1);
w_nb_n(3)=att_n_dot(3)+Cbn(3,1)*att_n_dot(1);

%%compute w_en_n & w_ie_n
if (styp==0)
    w_en_n=zeros(3,1);
    w_en_n(1)=vel_n(2)/(Re_eff);
    w_en_n(2)=-vel_n(1)/(Rn_eff);
    w_en_n(3)=-vel_n(2)*(sL/cL)/(Re_eff);

    w_ie_n=zeros(3,1);
    w_ie_n(1)=+WIE*cL;
    w_ie_n(2)=0;
    w_ie_n(3)=-WIE*sL;
elseif(styp==1)
    w_en_n=zeros(3,1);
    w_ie_n=zeros(3,1);
end

%attitude derivative
Cbn_dot=Cbn*skew(Cbn'*w_nb_n);

%%velocity derivative
vel_n_dot=Cbn*vel_b_dot+skew(w_nb_n)*vel_n;

%%position derivative
if (styp==0)
    pos_n_dot=zeros(3,1);
    pos_n_dot(1)=vel_n(1)/(Rn_eff);
    pos_n_dot(2)=vel_n(2)/(Re_eff)/cL;
    pos_n_dot(3)=-vel_n(3);
elseif (styp==1)
    pos_n_dot=zeros(3,1);
    pos_n_dot(1)=vel_n(1);
    pos_n_dot(2)=vel_n(2);
    pos_n_dot(3)=vel_n(3);
end

%%gyroscope output
gyro=Cbn'*(w_nb_n+w_en_n+w_ie_n);

% %%Acceleration output
% acc=Cbn'*(vel_ned_dot+(skew(w_en_n)+2*skew(w_ie_n))*vel_ned+gravity);

%%Acceleration output
vr_a=Cbn'*w_ie_n;
acc=vel_b_dot+skew(vr_a+gyro)*vel_b+Cbn'*gravity;

function [gyro_vib, vib_state_new, imu_vib_new]=update_vib(vib_state, vib_model, imu_vib, imu_or, dt)
nimu=size(vib_model,2);
vib_state_new=zeros(3,nimu);
imu_vib_new=zeros(3,3,nimu);
gyro_vib=zeros(3,nimu);
for in=1:nimu
    %compute vibration induced rotation rate in master frame
    Crig=imu_or(:,:,in);    %rigid Csm
    Cvib=imu_vib(:,:,in);   %vibration Cs's , complete rotation Cs'm=Crig*Cvib
    vib_state_new(:,in)=vib_model(1,in)*vib_state(:,in)+[0;0; vib_model(2,in)*randn(1)];
    gyro_vib(:,in)=Crig*(vib_state_new(:,in)-vib_state(:,in))/dt;

    %update the vibration
    rot=Cvib'*(vib_state_new(:,in)-vib_state(:,in));
    rot_norm=norm(rot);
    sr_a=1-(rot_norm^2/6)+(rot_norm^4/120);
    sr_b=(1/2)-(rot_norm^2/24)+(rot_norm^4/720);
    mx_a=eye(3)+sr_a*skew(rot)+sr_b*skew(rot)*skew(rot);
    Cvib=Cvib*mx_a;

    imu_vib_new(:,:,in)=Cvib;
end

repo: ekfmonoslam


text:
function y = cvx_value( x )
global cvx___
y = cellfun( @cvx_value, x, 'UniformOutput', false );

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [A,B,C,D,E] = bay_lssvm(model,level,type, nb, bay)
% Compute the posterior cost for the 3 levels in Bayesian inference
%
% >> cost = bay_lssvm({X,Y,type,gam,sig2}, level, type)
% >> cost = bay_lssvm(model              , level, type)
%
% Description
% Estimate the posterior probabilities of model (hyper-) parameters
% on the different inference levels:
%     - First level: In the first level one optimizes the support values alpha 's and the bias b.
%     - Second level: In the second level one optimizes the regularization parameter gam.
%     - Third level: In the third level one optimizes the kernel
%                    parameter. In the case of the common 'RBF_kernel' the kernel
%                    parameter is the bandwidth sig2.
%
% By taking the negative logarithm of the posterior and neglecting all constants, one
% obtains the corresponding cost. Computation is only feasible for
% one dimensional output regression and binary classification
% problems. Each level has its different in- and output syntax.
%
%
% Full syntax
%
%     1. Outputs on the first level
%
% >> [costL1,Ed,Ew,bay] = bay_lssvm({X,Y,type,gam,sig2,kernel,preprocess}, 1)
% >> [costL1,Ed,Ew,bay] = bay_lssvm(model, 1)
%
%       costL1 : Cost proportional to the posterior
%       Ed(*)  : Cost of the fitting error term
%       Ew(*)  : Cost of the regularization parameter
%       bay(*) : Object oriented representation of the results of the Bayesian inference
%
%     2. Outputs on the second level
%
% >> [costL2,DcostL2, optimal_cost, bay] = bay_lssvm({X,Y,type,gam,sig2,kernel,preprocess}, 2)
% >> [costL2,DcostL2, optimal_cost, bay] = bay_lssvm(model, 2)
%
%       costL2     : Cost proportional to the posterior on the second level
%       DcostL2(*) : Derivative of the cost
%       optimal_cost(*) : Optimality of the regularization parameter (optimal = 0)
%       bay(*)     : Object oriented representation of the results of the Bayesian inference
%
%     3. Outputs on the third level
%
% >> [costL3,bay] = bay_lssvm({X,Y,type,gam,sig2,kernel,preprocess}, 3)
% >> [costL3,bay] = bay_lssvm(model, 3)
%
%       costL3 : Cost proportional to the posterior on the third level
%       bay(*) : Object oriented representation of the results of the Bayesian inference
%
%     4. Inputs using the functional interface
%
% >> bay_lssvm({X,Y,type,gam,sig2,kernel,preprocess}, level)
% >> bay_lssvm({X,Y,type,gam,sig2,kernel,preprocess}, level, type)
% >> bay_lssvm({X,Y,type,gam,sig2,kernel,preprocess}, level, type, nb)
%
%         X            : N x d matrix with the inputs of the training data
%         Y            : N x 1 vector with the outputs of the training data
%         type         : 'function estimation' ('f') or 'classifier' ('c')
%         gam          : Regularization parameter
%         sig2         : Kernel parameter (bandwidth in the case of the 'RBF_kernel')
%         kernel(*)    : Kernel type (by default 'RBF_kernel')
%         preprocess(*) : 'preprocess'(*) or 'original'
%         level        : 1, 2, 3
%         type(*)      : 'svd'(*), 'eig', 'eigs', 'eign'
%         nb(*)        : Number of eigenvalues/eigenvectors used in the eigenvalue decomposition approximation
%
%     5. Inputs using the object oriented interface
%
% >> bay_lssvm(model, level, type, nb)
%
%         model    : Object oriented representation of the LS-SVM model
%         level    : 1, 2, 3
%         type(*)  : 'svd'(*), 'eig', 'eigs', 'eign'
%         nb(*)    : Number of eigenvalues/eigenvectors used in the eigenvalue decomposition approximation
%
%
% See also:
%   bay_lssvmARD, bay_optimize, bay_modoutClass, bay_errorbar

% Copyright (c) 2002,  KULeuven-ESAT-SCD, License & help @ http://www.esat.kuleuven.ac.be/sista/lssvmlab

%
% initiate and ev. preprocess
%
if ~isstruct(model), model = initlssvm(model{:}); end
model = prelssvm(model);
if model.y_dim>1,
  error(['Bayesian framework restricted to 1 dimensional regression' ...
	 ' and binary classification tasks']);
end

%
% train with the matlab routines
%model = adaptlssvm(model,'implementation','MATLAB');

eval('nb;','nb=ceil(sqrt(model.nb_data));');

if ~(level==1 | level==2 | level==3),
  error('level must be 1, 2 or 3.');
end

%
% delegate functions
%
if level==1,

  eval('type;','type=''train'';');
  %[cost, ED, EW, bay, model] = lssvm_bayL1(model, type);
  eval('[A,B,C,D,E] = lssvm_bayL1(model,type,nb,bay);','[A,B,C,D,E] = lssvm_bayL1(model,type,nb);');

elseif level==2,

  % default type
  eval('type;','type=''svd'';');
  %[costL2, DcostL2, optimal, bay, model] = lssvm_bayL2(model, type);

  eval('[A,B,C,D,E] = lssvm_bayL2(model,type,nb,bay);',...
       '[A,B,C,D,E] = lssvm_bayL2(model,type,nb);')

elseif level==3,

  % default type
  eval('type;','type=''svd'';');
  %[cost, bay, model] = lssvm_bayL3(model, bay);
  [A,B,C] = lssvm_bayL3(model,type,nb);

end

%
%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  FIRST LEVEL                   %
%                                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
function [cost, Ed, Ew, bay, model] = lssvm_bayL1(model, type, nb, bay)
%
% [Ed, Ew, cost,model] = lssvm_bayL1(model)
% [bay,model] = lssvm_bayL1(model)
%
% type = 'retrain', 'train', 'svd'
%
%

if ~(strcmpi(type,'train') | strcmpi(type,'retrain') | strcmpi(type,'eig') | strcmpi(type,'eigs')| strcmpi(type,'svd')| strcmpi(type,'eign')),
  error('type should be ''train'', ''retrain'', ''svd'', ''eigs'' or ''eign''.');
end
%type(1)=='t'
%type(1)=='n'

N = model.nb_data;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute Ed, Ew en costL1 based on training solution %
% TvG, Financial Timeseries Prediction using LS-SVM, 27-28 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (type(1)=='t'), % train
  % find solution of ls-svm
  model = trainlssvm(model);
  % prior %
  if model.type(1) == 'f',
    Ew = .5*sum(model.alpha.*  (model.ytrain(1:model.nb_data,:) - model.alpha./model.gam - model.b));
  elseif model.type(1) == 'c',
    Ew = .5*sum(model.alpha.*model.ytrain(1:model.nb_data,:).*  ...
		((1-model.alpha./model.gam)./model.ytrain(1:model.nb_data,:) - model.b));
  end

  % likelihood
  Ed = .5.*sum((model.alpha./model.gam).^2);

  % posterior
  cost = Ew+model.gam*Ed;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute Ed, Ew en costL1 based on SVD or nystrom %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
else

  if nargin<4,
    [bay.eigvals, bay.scores, ~, omega_r] = kpca(model.xtrain(model.selector,1:model.x_dim), ...
                                                  model.kernel_type, model.kernel_pars, [],type,nb,'original');

    bay.eigvals = bay.eigvals.*(N-1);
    bay.tol = 1000*eps;
    bay.Peff = find(bay.eigvals>bay.tol);
    bay.Neff = length(bay.Peff);
    bay.eigvals = bay.eigvals(bay.Peff);
    bay.scores = bay.scores(:,bay.Peff);
    %Zc = eye(N)-ones(model.nb_data)/model.nb_data;

    %disp('rescaling the scores');
    for i=1:bay.Neff,
      bay.Rscores(:,i) = bay.scores(:,i)./sqrt(bay.scores(:,i)'*bay.eigvals(i)*bay.scores(:,i));
    end
  end
  Y = model.ytrain(model.selector,1:model.y_dim);

  %%% Ew %%%%
  % (TvG: 4.75 - 5.73))
  YTM = (Y'-mean(Y))*bay.scores;
  Ew = .5*(YTM*diag(bay.eigvals)*diag((bay.eigvals+1./model.gam).^-2))*YTM';

  %%% cost %%%
  YTM = (Y'-mean(Y));
  %if model.type(1) == 'c', % 'classification'  (TvG: 5.74)
  %  cost = .5*YTM*[diag(bay.eigvals); zeros(model.nb_data-bay.Neff,bay.Neff)]*diag((bay.eigvals+1./model.gam).^-1)*bay.scores'*YTM';
  %elseif model.type(1) == 'f', % 'function estimation' % (TvG: 4.76)
			       % + correctie of zero eignwaardes
    cost = .5*(YTM*model.gam*YTM')-.5*YTM*bay.scores*diag((1+1./(model.gam.*bay.eigvals)).^-1*model.gam)*bay.scores'*YTM';
  %end

  %%% Ed %%%
  Ed = (cost-Ew)/model.gam;

end

bay.costL1 = cost;
bay.Ew = Ew;
bay.Ed = Ed;
bay.mu = (N-1)/(2*bay.costL1);
bay.zeta = model.gam*bay.mu;

% SECOND LEVEL
%
%
function [costL2, DcostL2, optimal, bay, model] = lssvm_bayL2(model,type,nb,bay)
%
%
%

if ~(strcmpi(type,'eig') | strcmpi(type,'eigs')| strcmpi(type,'svd')| strcmpi(type,'eign')),
  error('The used type needs to be ''svd'', ''eigs''  or ''eign''.')
end

  N = model.nb_data;
  % bayesian interference level 1

  eval('[cost, Ed, Ew, bay, model] = bay_lssvm(model,1,type,nb,bay); ',...
       '[cost, Ed, Ew, bay, model] = bay_lssvm(model,1,type,nb);');

  all_eigvals = zeros(N,1); all_eigvals(bay.Peff) = bay.eigvals;

  % Number of effective parameters
  bay.Geff = 1 + sum(model.gam.*all_eigvals ./(1+model.gam.*all_eigvals));
  bay.mu = .5*(bay.Geff-1)/(bay.Ew);
  bay.zeta = .5*(N-bay.Geff)/bay.Ed;
  % ideally: bay.zeta = model.gam*bay.mu;

  % log posterior (TvG: 4.73 - 5.71)
  costL2 = sum(log(all_eigvals+1./model.gam)) + (N-1).*log(bay.Ew+model.gam*bay.Ed);

  % gradient (TvG: 4.74 - 5.72)
  DcostL2 = -sum(1./(all_eigvals.*(model.gam.^2)+model.gam)) ...
	    + (N-1)*(bay.Ed/(bay.Ew+model.gam*bay.Ed));

  % endcondition fullfilled if optimal == 0;
  optimal = model.gam  - (N-bay.Geff)/(bay.Geff-1) * bay.Ew/bay.Ed;

  % update structure bay
  bay.optimal = optimal;
  bay.costL2 = costL2;
  bay.DcostL2 = DcostL2;

% THIRD LEVEL
%
%
function [costL3, bay, model] = lssvm_bayL3(model,type,nb)
%
% costL3 = lssvm_bayL3(model, type)
%

if ~(strcmpi(type,'svd') | strcmpi(type,'eigs') | strcmpi(type,'eign')),
  error('The used type needs to be ''svd'', ''eigs'' or ''eign''.')
end

% lower inference levels;
[model,~, bay] = bay_optimize(model,2,type,nb);

% test Neff << N
N = model.nb_data;
if sqrt(N)>bay.Neff,
  %model.kernel_pars
  %model.gam
  warning on;
  warning(['Number of degrees of freedom not tiny with respect to' ...
	   ' the number of datapoints. The approximation is not very good.']);
  warning off
end

% construct all eigenvalues
all_eigvals = zeros(N,1);
all_eigvals(bay.Peff) = bay.eigvals;

% L3 cost function
%costL3 = sqrt(bay.mu^bay.Neff*bay.zeta^(N-1)./((bay.Geff-1)*(N-bay.Geff)*prod(bay.mu+bay.zeta.*all_eigvals)));
%costL3 = .5*bay.costL2 - log(sqrt(2/(bay.Geff-1))) - log(sqrt(2/(N-bay.Geff)))
costL3 = -(bay.Neff*log(bay.mu) + (N-1)*log(bay.zeta)...
	 - log(bay.Geff-1) -log(N-bay.Geff) - sum(log(bay.mu+bay.zeta.*all_eigvals)));
bay.costL3 = costL3;

author: peterhcharlton

repo: RRest


text:
function opts = U_param(opts,mn,ln)

opts.do_thres=1; % threshold the output image [0 255] after each iteration
opts.grad_ran=[0.1 100]; % adaptively scale the gradient scale if it's out of range
switch opts.task
    case 0 % feature inversion
        switch mn
            case 'alex'
            opts.lambdaL2 = 8e-10 ;opts.lambdaTV = 1e0;opts.lambdaD = 1e2;
            opts.learningRate = 0.0004 * [...
                ones(1,100), ...
                0.1 * ones(1,200), ...
                0.01 * ones(1,100)];
              switch ln
                case 'p5'; % default
                opts.lambdaL2 = 8e-8 ;opts.lambdaTV = 1e2;opts.lambdaD = 1e1;
            opts.learningRate = 0.0004 * [...
                ones(1,500), ...
                0.1 * ones(1,200), ...
                0.01 * ones(1,100)];
                end
            end
    case 1 % neuron inversion
        switch mn
            case 'alex'
                opts.lambdaL2 = 1e-10 ;opts.lambdaTV = 1e2;opts.lambdaD = 1e-3;
                opts.learningRate = 0.04 * [...
                    0.05 * ones(1,100), ...
                    0.01 * ones(1,100)];
                switch ln
                case 'p5'; % default
                case 'f7'; opts.lambdaD = 3e-2;
                case 'f8'; opts.lambdaD = 1e-2;
                end
            case 'vgg16'
                opts.lambdaL2 = 1e-10 ;opts.lambdaTV = 8e2;opts.lambdaD = 3e-3;
                opts.learningRate = 0.001 * [...
                    0.05 * ones(1,100), ...
                    0.01 * ones(1,50)];
                opts.grad_ran=[0.1 100];
            case 'nin'
                opts.lambdaL2 = 1e-10 ;opts.lambdaTV = 8e2;opts.lambdaD = 3e-3;
                opts.learningRate = 0.004 * [...
                    0.05 * ones(1,100), ...
                    0.01 * ones(1,50)];
            case 'gnet'
                % annealing schedule
                opts.lambdaL2 = [1e0*ones(1,50) 1e1*ones(1,50) 1e2*ones(1,100)];
                opts.lambdaTV = 1e-1*[5e2*ones(1,100) 5e2*ones(1,100)];
                opts.lambdaD = 3e-3;
                opts.learningRate = 0.04 * [...
                    0.05 * ones(1,100), ...
                    0.01 * ones(1,100)];
        end
    case 2
        % image inpainting
        switch mn
        case 'alex'
            opts.lambdaL2 = 1e-8 ;opts.lambdaTV = 5e2;
            opts.lambdaD = 5e1;
            opts.grad_ran=[0 20];
            opts.learningRate = 0.004 * [...
            0.05 * ones(1,150), ...
            0.05 * ones(1,0), ...
            0.05 * ones(1,0)];
        end
end

author: donglaiw

repo: mNeuron

save_path: github-repos/MATLAB/donglaiw-mNeuron


text:
function [fx] = evolution0bisND(x,theta,u,in)
% 0-ToM's evolution function (without doubled hidden states)
% function [fx] = evolution0bisND(x,theta,u,in)
% 0-ToM is simply tracking the log-odds of P(o=1), where o is the
% opponent's action. This variable is updated according to a Laplace-Kalman
% filter, yielding 2 sufficient statistics, m and V. In this scheme, the
% only evolution param (theta) is 0-ToM's prior volatity about her
% opponent's log-odds.
% IN:
%   - x: sufficient statistics of log-odds of P(o=1):
%       x(1)= E[log-odds]
%       x(2)= log V[log-odds] (log-scale for numerical reasons)
%   - theta: 0-ToM's prior (log-) volatity
%   - u: u(1)= last opponent's move (o)
%   - in: [useless here]
% OUT:
%   - fx: updated sufficient statistics of log-odds of P(o=1)

if isempty(u)||isnan(u(1)) % missed trial
    fx = x; % no update
else % trial OK
    % -- deal with superceding competition with other generative models --
    % [See, e.g., f_metaTom.m]
    try
        w = inF.metaweight;
    catch
        w = 1;
    end
    % -- learning rule --
    m0 = x(1); % current E[log-odds]
    V0 = exp(x(2)); % current V[log-odds]
    p0 = VBA_sigmoid(m0); % current estimate of P(o=1)
    volatility = exp(theta(1));
    V = 1./((1./(volatility+V0))+w*p0*(1-p0)); % updated V[log-odds]
    m = m0 + w*V*(u(1)-p0); % updated E[log-odds] (Laplace-Kalman update rule)
    % wrap-up
    fx = [VBA_sigmoid(VBA_sigmoid(m),'inverse',true);log(V)]; % for numerical purposes
end

author: MBB-team

repo: VBA-toolbox


text:
function auc = aucroc(p_predicted, p_target, freq)

        % Count observations by class
        nTarget     = sum(freq .* p_target);
        nBackground = sum(freq .* (1-p_target));

        % Rank data
        R = tiedrank(p_predicted);  % 'tiedrank' from Statistics Toolbox
        [R_uniq, I, J] = unique(R);
        [R_mean,R_num] = grpstats(freq,R,{'mean','numel'});
        R_freq = R_mean .* R_num;

        % adjust for counts
        rank_start_freq = cumsum(R_freq);
        rank_start_freq = [0; rank_start_freq(1:(end-1))] + 1;
        rank_mean_freq = (2*rank_start_freq+R_freq-1)/2;

        rank_orig = rank_mean_freq(J);

        % Calculate AUC
        %Error = (sum(R(Actual == 1)) - (nTarget^2 + nTarget)/2) / (nTarget * nBackground);
        auc = (sum(rank_orig .* p_target .* freq) - (nTarget^2 + nTarget)/2) / (nTarget * nBackground);

author: Sable

repo: mcbench-benchmarks


text:
function t = subsasgn(t,s,b)
%SUBSASGN Subscripted assignement for ktensor.
%
%   See also KTENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

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

switch s(1).type
    case '.'
        switch s(1).subs
            case 'lambda'
                if length(s) == 1
                    t = ktensor(b, t.u);
                else
                    newlambda = subsasgn(t.lambda, s(2:end), b);
                    t = ktensor(newlambda, t.u);
                end
            case {'u','U'}
                if length(s) == 1
                    t = ktensor(t.lambda, b);
                else
                    tmpu = subsasgn(t.u, s(2:end), b);
                    t = ktensor(t.lambda, tmpu);
                end
            otherwise
                error(['No such field: ', s(1).subs]);
        end
    case '()'
        error('Cannot change individual entries in a ktensor.')
    case '{}'
        new_s(1).type = '.';
        new_s(1).subs = 'u';
        new_s(2:length(s)+1) = s;
        t = subsasgn(t, new_s, b);
    otherwise
        error('Invalid subsasgn.');
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [interpV] = finterp3(xInterpV, yInterpV, zInterpV, field3M, xFieldV, yFieldV, zFieldV, OOBV)
%"finterp3"
%   Fast 3D linear interpolation of the data field3M defined at
%   xFieldV,yFieldV,zFieldV, to the points listed in x/y/zInterpV.
%
%   xFieldV and yFieldV are 3 element vectors of the form [start, delta, end],
%   where start is the coordinate of the first column(x) or row(y) and end
%   is the coord of the last column/row, with delta being the spacing.  Delta
%   can be negative if start > end.
%
%   zFieldV must be specified as a list of coordinates for each slice in
%   field3M.
%
%   xInterpV, yInterpV, zInterpV are vectors of equal length indicating the
%   x,y,z, coordinates to interpolate from original data field3M.
%
%   Based on code by J.O.Deasy.
%
%JRA 1/13/05
%
%Usage:
%   function [interpV] = finterp3(xInterpV, yInterpV, zInterpV, field3M, xFieldV, yFieldV, zFieldV, outOfBoundsVal)

siz = size(field3M);

if ~exist('OOBV')
    OOBV = NaN;
end

xDelta = xFieldV(2);
yDelta = yFieldV(2);

%Get r,c,s indices.
% cols = interp1q(reshape(xFieldV, [], 1), reshape(1:length(xFieldV),[], 1), reshape(xInterpV, [], 1));
% rows = interp1q(reshape(yFieldV, [], 1), reshape(1:length(yFieldV),[], 1), reshape(yInterpV, [], 1));
cols = (xInterpV-(xFieldV(1)-xDelta))/xDelta;
rows = (yInterpV-(yFieldV(1)-yDelta))/yDelta;
if length(zFieldV) > 1;
    slcs = interp1(zFieldV,1:length(zFieldV),zInterpV);
else
    slcs = ones(size(cols));    %This effectively negates Z.  All values are in plane.  Bad idea?
end

clear xInterpV yInterpV zInterpV
clear xFieldV yFieldV zFieldV

%Find indices out of bounds.
colNaN = cols >= siz(2) | cols < 1;
rowNaN = rows >= siz(1) | rows < 1;
% slcNaN = slcs >= siz(3) | slcs < 1;

% colNaN = isnan(cols);
% rowNaN = isnan(rows);
slcNaN = isnan(slcs) | slcs < 1 | slcs >= siz(3);

%Set those to a proxy 1.
rows(rowNaN) = 1;
cols(colNaN) = 1;
slcs(slcNaN) = 1;

colFloor = floor(cols);
colMod   = cols - colFloor;
oneMinusColMod = (1-colMod);

rowFloor = floor(rows);
rowMod = rows - rowFloor;
oneMinusRowMod = (1-rowMod);

slcFloor = floor(slcs);
slcMod = slcs - slcFloor;
oneMinusSlcMod = (1-slcMod);

clear rows cols slcs

%Linear indices of lower bound contributing points.
INDEXLIST = rowFloor + (colFloor-1)*siz(1) + (slcFloor-1)*siz(1)*siz(2);

clear rowFloor colFloor slcFloor

%Linear offsets when moving in 3D matrix.
oneRow = 1;
oneCol = siz(1);
oneSlc = siz(1)*siz(2);

%Accumulate contribution from each voxel surrounding x,y,z point.
interpV = double(field3M(INDEXLIST)) .* oneMinusRowMod .* oneMinusColMod .* oneMinusSlcMod;
interpV = interpV + double(field3M(INDEXLIST+(oneRow))) .* rowMod .* oneMinusColMod .* oneMinusSlcMod;
interpV = interpV + double(field3M(INDEXLIST+(oneCol))) .* oneMinusRowMod .* colMod .* oneMinusSlcMod;
interpV = interpV + double(field3M(INDEXLIST+(oneCol+oneRow))) .* rowMod .* colMod .* oneMinusSlcMod;;
interpV = interpV + double(field3M(INDEXLIST+(oneSlc))) .* oneMinusRowMod .* oneMinusColMod .* slcMod;
interpV = interpV + double(field3M(INDEXLIST+(oneSlc+oneRow))) .* rowMod .* oneMinusColMod .* slcMod;
interpV = interpV + double(field3M(INDEXLIST+(oneSlc+oneCol))) .* oneMinusRowMod .* colMod .* slcMod;
interpV = interpV + double(field3M(INDEXLIST+(oneSlc+oneCol+oneRow))) .* rowMod .* colMod .* slcMod;

%Replace proxy 1s with out of bounds vals.
interpV(rowNaN | colNaN | slcNaN) = OOBV;

return;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [hess,pattern,symhess] = symHess(fun,nvar,isTril,var)
%SYMHESS  Returns a Symbolically Differentiated Hessian
%
%   hess = symHess(fun) uses the symbolic toolbox to automatically generate
%   the Hessian of the function handle fun. You should supply the function
%   handle in the form @(x) x(1)^2 + x(2), noting x must be the only
%   variable used, and is indexed in the equation (no vector operations).
%
%   hess = symHess(fun,nvar) specifies the number of variables in the
%   equation, assuming consecutive ordering. Useful if a variable is not
%   specified in the original equation to pad the Hessian with zeros.
%
%   hess = symHess(fun,nvar,isTril) specifies if the returned Hessian
%   should be Symmetric Lower Triangular.
%
%   [hess,pattern] = symHess(...) also returns the sparsity pattern as a
%   function handle.

%   Copyright (C) 2013 Jonathan Currie (I2C2)

if(nargin < 4), var = 'x'; end
if(nargin < 3), isTril = false; end
if(nargin < 2 || isempty(nvar)), nvar = 0; end

if(~optiCheckSymTBX())
    hess = []; symhess = [];
    return
end

if(~isa(fun,'function_handle') && ~isa(fun,'barvec'))
    error('Fun should be a function handle!');
end
if(isa(fun,'function_handle') && nargin(fun) ~= 1)
    error('fun should only have one input argument (x)');
end

%Convert function to a symbolic expression
if(isa(fun,'function_handle'))
    [symfun,ind] = func2sym(fun,var);
    ind = ind{1};
else
    symfun = sym(getEq(fun)); %convert from barvec to symbolic expression
    v = char(symvar(symfun));
    ind = 1:length(strfind(v,var));
end

%Check Symbolic Expression is Scalar
if(~isscalar(symfun))
    error('symHess only works on scalar functions');
end

%Create each symbolic partial derivative
symhess = symPartialDer(symfun,var,nvar,ind,true);

%Convert if required to lower triangular
if(isTril), symhess = tril(symhess); end

%Extract Sparsity Pattern
sp = zeros(size(symhess)); sp(symhess ~= 0) = 1; sp = sparse(sp);
pattern = @() sp;

%Return to function handle
hess = sym2func(symhess,var);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [pvec, pstruct] = tapas_hgf_binary_mab_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/>.

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

l = r.c_prc.n_levels;

pvec(1:l)         = ptrans(1:l);                                 % mu_0
pstruct.mu_0      = pvec(1:l);
pvec(l+1:2*l)     = exp(ptrans(l+1:2*l));                        % sa_0
pstruct.sa_0      = pvec(l+1:2*l);
pvec(2*l+1:3*l)   = ptrans(2*l+1:3*l);                           % rho
pstruct.rho       = pvec(2*l+1:3*l);
pvec(3*l+1:4*l-1) = exp(ptrans(3*l+1:4*l-1));                    % ka
pstruct.ka        = pvec(3*l+1:4*l-1);
pvec(4*l:5*l-1)   = ptrans(4*l:5*l-1);                           % om
pstruct.om        = pvec(4*l:5*l-1);

return;

author: translationalneuromodeling

repo: tapas


text:
function prob = yalmip2mosek(interfacedata);

% Retrieve needed data
options = interfacedata.options;
F_struc = interfacedata.F_struc;
c       = interfacedata.c;
Q       = interfacedata.Q;
K       = interfacedata.K;
x0      = interfacedata.x0;
integer_variables = interfacedata.integer_variables;
binary_variables = interfacedata.binary_variables;
extended_variables = interfacedata.extended_variables;
ub      = interfacedata.ub;
lb      = interfacedata.lb;
mt      = interfacedata.monomtable;

% *********************************
% What type of variables do we have
% *********************************
linear_variables = find((sum(abs(mt),2)==1) & (any(mt==1,2)));
nonlinear_variables = setdiff((1:size(mt,1))',linear_variables);
sigmonial_variables = find(any(0>mt,2) | any(mt-fix(mt),2));

if ~isempty(sigmonial_variables) | isequal(interfacedata.solver.version,'GEOMETRIC')
    prob = create_mosek_geometric(options,F_struc,c,Q,K,ub,lb,mt,linear_variables,extended_variables);
else
    prob = create_mosek_lpqp(options,F_struc,c,Q,K,ub,lb,mt,linear_variables,integer_variables);
end

function prob = create_mosek_lpqp(options,F_struc,c,Q,K,ub,lb,mt,linear_variables,integer_variables);

prob.c = c;
if ~isempty(F_struc)
    prob.a = -F_struc(:,2:end);
    prob.buc = full(F_struc(:,1));
    prob.blc = repmat(-inf,length(prob.buc),1);
else
    prob.a = sparse(ones(1,length(c))); % Dummy constraint
    prob.buc = inf;
    prob.blc = -inf;
end
if isempty(lb)
    prob.blx = repmat(-inf,1,length(c));
else
    prob.blx = lb;
end
if isempty(ub)
    prob.bux = repmat(inf,1,length(c));
else
    prob.bux = ub;
end

if K.f>0
    prob.blc(1:K.f) = prob.buc(1:K.f);
end

[prob.qosubi,prob.qosubj,prob.qoval] = find(tril(sparse(2*Q)));

if any(K.q)
    nof_new = sum(K.q);
    prob.a = [prob.a [spalloc(K.f,nof_new,0);spalloc(K.l,nof_new,0);speye(nof_new)]];
    %prob.a(1+K.f+K.l:end,1:length(c)) = prob.a(1+K.f+K.l:end,1:length(c));
    prob.blc(1+K.f+K.l:end) = prob.buc(1+K.f+K.l:end);
    prob.buc(1+K.f+K.l:end) = prob.buc(1+K.f+K.l:end);
    prob.c = [prob.c;zeros(nof_new,1)];
    top = size(F_struc,2)-1;
    for i = 1:length(K.q)
        prob.cones{i}.type = 'MSK_CT_QUAD';
        prob.cones{i}.sub  = top+1:top+K.q(i);
        prob.blx(top+1:top+K.q(i)) = -inf;
        prob.bux(top+1:top+K.q(i)) = inf;
        top = top + K.q(i);
    end
end

if ~isempty(integer_variables)
    prob.ints.sub = integer_variables;
end

prob.param = options.mosek;

function prob = create_mosek_geometric(options,F_struc,c,Q,K,ub,lb,mt,linear_variables,extended_variables);

x = [];
D_struc = [];
res = [];
solvertime = 0;
[prob,problem] = yalmip2geometric(options,F_struc,c,Q,K,ub,lb,mt,linear_variables,extended_variables);
if problem == 0

    % Mosek does not support equalities
    if ~isempty(prob.G)
        prob.A = [prob.A;prob.G;-prob.G];
        prob.b = [prob.b;prob.h;1./prob.h];
        prob.map = [prob.map;max(prob.map) + (1:2*length(prob.h))'];
    end

else
    prob = [];
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function net = init_LapSRN_model(opts)
% -------------------------------------------------------------------------
%   Description:
%       initialize LapSRN model
%
%   Input:
%       - opts  : options generated from init_LapSRN_opts()
%
%   Output:
%       - net   : dagnn model
%
%   Citation:
%       Deep Laplacian Pyramid Networks for Fast and Accurate Super-Resolution
%       Wei-Sheng Lai, Jia-Bin Huang, Narendra Ahuja, and Ming-Hsuan Yang
%       IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2017
%
%   Contact:
%       Wei-Sheng Lai
%       wlai24@ucmerced.edu
%       University of California, Merced
% -------------------------------------------------------------------------

    %% parameters
    rng('default');
    rng(0) ;

    f       = opts.conv_f;
    n       = opts.conv_n;
    pad     = floor(f/2);
    depth   = opts.depth;
    scale   = opts.scale;
    level   = ceil(log(scale) / log(2));
    if( f == 3 )
        crop = [0, 1, 0, 1];
    elseif( f == 5 )
        crop = [1, 2, 1, 2];
    else
        error('Need to specify crop in deconvolution for f = %d\n', f);
    end

    %% initialize model
    net = dagnn.DagNN;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Feature extraction branch
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    sigma   = opts.init_sigma;
    filters = sigma * randn(f, f, 1, n, 'single');
    biases  = zeros(1, n, 'single');

    % conv
    inputs  = { 'LR' };
    outputs = { 'input_conv' };
    params  = { 'input_conv_f', 'input_conv_b' };

    net.addLayer(outputs{1}, ...
                 dagnn.Conv('size', size(filters), ...
                            'pad', pad, ...
                            'stride', 1), ...
                 inputs, outputs, params);

    idx = net.getParamIndex(params{1});
    net.params(idx).value         = filters;
    net.params(idx).learningRate  = 1;
    net.params(idx).weightDecay   = 1;

    idx = net.getParamIndex(params{2});
    net.params(idx).value         = biases;
    net.params(idx).learningRate  = 0.1;
    net.params(idx).weightDecay   = 1;

    % ReLU
    inputs  = { 'input_conv' };
    outputs = { 'input_relu' };

    net.addLayer(outputs{1}, ...
                 dagnn.ReLU('leak', 0.2), ...
                 inputs, outputs);

    next_input = outputs{1};

    %% deep conv layers (f x f x n x n)
    sigma   = sqrt( 2 / (f * f * n) );

    for s = level : -1 : 1

        % conv layers (f x f x n x n)
        for d = 1:depth

            filters = sigma * randn(f, f, n, n, 'single');
            biases  = zeros(1, n, 'single');

            % conv
            inputs  = { next_input };
            outputs = { sprintf('level%d_conv%d', s, d) };
            params  = { sprintf('level%d_conv%d_f', s, d), ...
                        sprintf('level%d_conv%d_b', s, d)};

            net.addLayer(outputs{1}, ...
                         dagnn.Conv('size', size(filters), ...
                                    'pad', pad, ...
                                    'stride', 1), ...
                         inputs, outputs, params);

            idx = net.getParamIndex(params{1});
            net.params(idx).value         = filters;
            net.params(idx).learningRate  = 1;
            net.params(idx).weightDecay   = 1;

            idx = net.getParamIndex(params{2});
            net.params(idx).value         = biases;
            net.params(idx).learningRate  = 0.1;
            net.params(idx).weightDecay   = 1;

            % ReLU
            inputs  = { sprintf('level%d_conv%d', s, d) };
            outputs = { sprintf('level%d_relu%d', s, d) };

            net.addLayer(outputs{1}, ...
                         dagnn.ReLU('leak', 0.2), ...
                     inputs, outputs);

            next_input = outputs{1};

        end

        %% features upsample layers
        filters = sigma * randn(f, f, n, n, 'single');
        biases  = zeros(1, n, 'single');

        inputs  = { next_input };
        outputs = { sprintf('level%d_upconv', s) };
        params  = { sprintf('level%d_upconv_f', s), ...
                    sprintf('level%d_upconv_b', s) };

        net.addLayer(outputs{1}, ...
                     dagnn.ConvTranspose(...
                         'size', size(filters), ...
                         'upsample', 2, ...
                         'crop', crop, ...
                         'numGroups', 1, ...
                         'hasBias', true), ...
                     inputs, outputs, params) ;

        idx = net.getParamIndex(params{1});
        net.params(idx).value         = filters;
        net.params(idx).learningRate  = 1;
        net.params(idx).weightDecay   = 1;

        idx = net.getParamIndex(params{2});
        net.params(idx).value         = biases;
        net.params(idx).learningRate  = 0.1;
        net.params(idx).weightDecay   = 1;

        %% ReLU
        inputs  = { sprintf('level%d_upconv', s) };
        outputs = { sprintf('level%d_uprelu', s) };

        net.addLayer(outputs{1}, ...
                     dagnn.ReLU('leak', 0.2), ...
                     inputs, outputs);

        next_input = outputs{1};

        %% residual prediction layer (f x f x n x 1)
        sigma   = sqrt(2 / (f * f * n));
        filters = sigma * randn(f, f, n, 1, 'single');
        biases  = zeros(1, 1, 'single');

        inputs  = { next_input };
        outputs = { sprintf('level%d_residual', s) };
        params  = { sprintf('level%d_residual_conv_f', s), ...
                    sprintf('level%d_residual_conv_b', s) };

        net.addLayer(outputs{1}, ...
            dagnn.Conv('size', size(filters), ...
                       'pad', pad, ...
                       'stride', 1), ...
            inputs, outputs, params);

        idx = net.getParamIndex(params{1});
        net.params(idx).value         = filters;
        net.params(idx).learningRate  = 1;
        net.params(idx).weightDecay   = 1;

        idx = net.getParamIndex(params{2});
        net.params(idx).value         = biases;
        net.params(idx).learningRate  = 0.1;
        net.params(idx).weightDecay   = 1;

    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Image reconstruction branch
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    next_input = 'LR';

    for s = level : -1 : 1

        %% image upsample layer
        filters = single(bilinear_kernel(4, 1, 1));

        inputs  = { next_input };
        outputs = { sprintf('level%d_img_up', s) };
        params  = { sprintf('level%d_img_up_f', s) };

        net.addLayer(outputs{1}, ...
            dagnn.ConvTranspose(...
                'size', size(filters), ...
                'upsample', 2, ...
                'crop', 1, ...
                'numGroups', 1, ...
                'hasBias', false), ...
            inputs, outputs, params) ;

        idx = net.getParamIndex(params{1});
        net.params(idx).value         = filters;
        net.params(idx).learningRate  = 1;
        net.params(idx).weightDecay   = 1;

        %% residual addition layer
        inputs  = { sprintf('level%d_img_up', s), ...
                    sprintf('level%d_residual', s) };
        outputs = { sprintf('level%d_output', s) };
        net.addLayer(outputs{1}, ...
            dagnn.Sum(), ...
            inputs, outputs);

        next_input = outputs{1};

        %% Loss layer
        inputs  = { next_input, ...
                    sprintf('level%d_HR', s) };
        outputs = { sprintf('level%d_%s_loss', s, opts.loss) };

        net.addLayer(outputs{1}, ...
                 dagnn.vllab_dag_loss(...
                    'loss_type', opts.loss), ...
                 inputs, outputs);

    end

end

repo: LapSRN


text:
function x_opt = plotInternalModel(internalmodel,x,n,localindex,color,opts)
% Code used by both lmi/plot and optimizer/plot

if isempty(internalmodel.binary_variables)
    [x_opt{1},errorstatus] = generateBoundary(internalmodel,x,n,localindex);
else
    if strcmp(internalmodel.solver.tag,'BNB')
        internalmodel.solver = internalmodel.solver.lower;
    end
    nBin = length(internalmodel.binary_variables);
    p = extractLP(internalmodel);
    p = extractOnly(p,internalmodel.binary_variables);

    internalmodel.F_struc = [zeros(nBin,size(internalmodel.F_struc,2));internalmodel.F_struc];
    I = eye(nBin);
    internalmodel.F_struc(1:nBin,1+internalmodel.binary_variables) = I;
    internalmodel.K.f = internalmodel.K.f + length(internalmodel.binary_variables);
    errorstatus = 1;
    x_opt = {};
    errorstatus = zeros(1,2^nBin);
    for i = 0:2^nBin-1;
        comb = dec2decbin(i,nBin);
        if checkfeasiblefast(p,comb(:),1e-6)
            internalmodel.F_struc(1:nBin,1) = -comb(:);
            [x_temp,wrong] = generateBoundary(internalmodel,x,n,localindex);
            if ~wrong
                errorstatus(i+1) = 0;
                x_opt{end+1} = x_temp;
            end
        else
            errorstatus(i+1)=0;
        end
    end
end

if all(errorstatus)
    if nargout==0
        plot(0);
    end
elseif nargout == 0
    for i = 1:length(x_opt)
        try
            plotSet(x_opt{i},color(1+rem(i-1,size(color,1)),:),opts);
        catch
        end
    end
end

if nargout > 0
    varargout{1} = x_opt;
end

function [xout,errorstatus] = solvefordirection(c,internalmodel,uv)
internalmodel.c = 0*internalmodel.c;
internalmodel.c(uv) = c;
sol  = feval(internalmodel.solver.call,internalmodel);
xout = sol.Primal;
xout = xout(uv(:));
errorstatus = sol.problem;

function p = plotSet(x_opt,color,options)
if size(x_opt,1)==1
    p = line(x_opt,[0 0],'color',color);
    set(p,'LineStyle',options.plot.wirestyle);
    set(p,'LineStyle',options.plot.wirestyle);
    set(p,'LineWidth',options.plot.linewidth);
    set(p,'EdgeColor',options.plot.edgecolor);
    set(p,'Facealpha',options.plot.shade);
elseif size(x_opt,1)==2
    p = patch(x_opt(1,:),x_opt(2,:),color);
    set(p,'LineStyle',options.plot.wirestyle);
    set(p,'LineWidth',options.plot.linewidth);
    set(p,'EdgeColor',options.plot.edgecolor);
    set(p,'Facealpha',options.plot.shade);
else
    try
        K = convhulln(x_opt');
        p = patch('Vertices', x_opt','Faces',K,'FaceColor', color);
    catch
         p = fill3(x_opt(1,:),x_opt(2,:),x_opt(3,:),1);
    end
    set(p,'LineStyle',options.plot.wirestyle);
    set(p,'LineWidth',options.plot.linewidth);
    set(p,'EdgeColor',options.plot.edgecolor);
    set(p,'Facealpha',options.plot.shade);
    lighting gouraud;
    view(3);
    camlight('headlight','infinite');
    camlight('headlight','infinite');
    camlight('right','local');
    camlight('left','local');
end

function [x_opt,errorstatus] = generateBoundary(internalmodel,x,n,localindex);

x_opt = [];
phi = [];
errorstatus = 0;
waitbar_created = 0;
t0 = clock;
waitbar_starts_at = 2;
lastdraw = clock;
try % Try to ensure that we close h
    if length(x)==2
        mu = 0.5;
    else
        mu=1;
    end
    n_ = n;
    n = ceil(mu*n);
   % h = waitbar(0,['Please wait, solving ' num2str(n_) ' problems using ' internalmodel.solver.tag]);
    angles = (0:(n))*2*pi/n;
    if length(x)==2
        c = [cos(angles);sin(angles)];
    elseif length(x) == 1
        c = [-1 1];n = 2;
    else
        c = randn(3,n);
    end
    i=1;
    while i<=n & errorstatus ~=1
        [xi,errorstatus] = solvefordirection(c(:,i),internalmodel,localindex(:));
        if errorstatus == 2
            disp('Discovered unbounded direction. You should add bounds on variables')
        elseif errorstatus == 12
            [xi,errorstatus] = solvefordirection(0*c(:,i),internalmodel,localindex(:));
            if errorstatus == 0
                errorstatus = 2;
                disp('Discovered unbounded direction. You should add bounds on variables')
            end
        end
        x_opt = [x_opt xi];
        if ~waitbar_created
            if etime(clock,t0)>waitbar_starts_at;
             h = waitbar(0,['Please wait, solving ' num2str(n_) ' problems using ' internalmodel.solver.tag]);
             waitbar_created = 1;
            end
        end
        if waitbar_created & etime(clock,lastdraw)>1/10
            waitbar(i/n_,h)
            lastdraw = clock;
        end
        i=i+1;
    end

    if errorstatus==0 & length(x)==2
        % Close the set
        x_opt = [x_opt x_opt(:,1)];

        % Add points adaptively
        pick = 1;
        n = floor((1-mu)*n_);
        for i = 1:1:n
            for j= 1:(size(x_opt,2)-1)
                d = x_opt(:,j)-x_opt(:,j+1);
                distance(j,1) = d'*d;
            end
            [dist,pos]=sort(-distance);
            % Select insertion point
            phii=(angles(pos(pick))+angles(pos(pick)+1))/2;
            xi = solvefordirection([cos(phii);sin(phii)],internalmodel,localindex);
            d1=xi-x_opt(:,pos(pick));
            d2=xi-x_opt(:,pos(pick)+1);
            if d1'*d1<1e-3 | d2'*d2<1e-3
                pick = pick+1;
            else
                angles = [angles(1:pos(pick)) phii  angles((pos(pick))+1:end)];
                x_opt = [x_opt(:,1:pos(pick)) xi  x_opt(:,(pos(pick))+1:end)];
            end
            if ~waitbar_created
            if etime(clock,t0)>waitbar_starts_at;
                h = waitbar(0,['Please wait, solving ' num2str(n_) ' problems using ' internalmodel.solver.tag]);
                waitbar_created = 1;
            end
            end
            if waitbar_created
            waitbar((ceil(n_*mu)+i)/n_,h);
            end
        end
    end
    if waitbar_created
    close(h);
    end
catch
    if waitbar_created
    close(h);
    end
end

function pLP = extractLP(p)
pLP = p;
pLP.F_struc = pLP.F_struc(1:p.K.f+p.K.l,:);
pLP.K.q = 0;
pLP.K.s = 0;

function pRed = extractOnly(p,these)
pRed = p;
p.F_struc(:,1+these) = 0;
removeEQ = find(any(p.F_struc(1:pRed.K.f,2:end),2));
removeLP = find(any(p.F_struc(1+pRed.K.f:end,2:end),2));
pRed.F_struc(pRed.K.f+removeLP,:)=[];
pRed.F_struc(removeEQ,:)=[];
pRed.K.f = pRed.K.f - length(removeEQ);
pRed.K.l = pRed.K.l - length(removeLP);
pRed.F_struc = pRed.F_struc(:,[1 1+these]);
pRed.lb = pRed.lb(these);
pRed.ub = pRed.ub(these);


text:
function [w,A] = train_blr(X,y,var,Sigma_p)
%TRAIN_BLR Train Bayesian Linear Regression
%
%   Bayesian Linear Regression, we seek to learn a regression function
%   where the relation between the predictor variable, y and the domain X
%   is linear in the parameters whilst assuming a prior probability
%   distribtion over the parameters
%
%     y = f(x) + \epsilon, \epsilon ~ N(0,var) (regression model)    (1)
%
%     y = x' * w with w ~ N(0,\Sigma_p)         (2)
%
%   input ----------------------------------------------------------------
%
%       o X: (N x D), N samples of dimension D
%
%       o y: (D x 1), target value
%
%       o Sigma_p: (D x D), Prior Covariance matrix on the weights
%
%   output ----------------------------------------------------------------
%
%       o w: (D x 1), weight vectors of equation (1)-(2)
%
%
%
N = size(X,1);

%                               (D x 1)
%                        (D x D)        (D x N)  * (N x 1)
%             ((D x N) * (N x D) + (D x D))
X  = [X,ones(N,1)]; % add bias

A = (1/var) .* X'*X + inv(Sigma_p);
w = (1/var) .* (A \ (X' * y));

end

author: epfl-lasa


text:
% ********************************************************************
% * mex File compiling code for Random Forest (for linux)
% * mex interface to Andy Liaw et al.'s C code (used in R package randomForest)
% * Added by Abhishek Jaiantilal ( abhishek.jaiantilal@colorado.edu )
% * License: GPLv2
% * Version: 0.02
% ********************************************************************/

function compile_windows
    system('del *.mexw32;del *.mexw64;');

    fprintf('I am going to use the precompiled fortran file\n');
    fprintf('If it doesnt work then use cygwin+g77 (or gfortran) to recompile rfsub.f\n');

    if strcmp(computer,'PCWIN64')
        mex  -DMATLAB -DWIN64 -output mexClassRF_train   src/classRF.cpp src/classTree.cpp src/cokus.cpp precompiled_rfsub/win64/rfsub.o src/mex_ClassificationRF_train.cpp   src/rfutils.cpp
        mex  -DMATLAB -DWIN64 -output mexClassRF_predict src/classRF.cpp src/classTree.cpp src/cokus.cpp precompiled_rfsub/win64/rfsub.o src/mex_ClassificationRF_predict.cpp src/rfutils.cpp
    elseif strcmp(computer,'PCWIN')
        mex  -DMATLAB -output mexClassRF_train   src/classRF.cpp src/classTree.cpp src/cokus.cpp precompiled_rfsub/win32/rfsub.o src/mex_ClassificationRF_train.cpp   src/rfutils.cpp
        mex  -DMATLAB -output mexClassRF_predict src/classRF.cpp src/classTree.cpp src/cokus.cpp precompiled_rfsub/win32/rfsub.o src/mex_ClassificationRF_predict.cpp src/rfutils.cpp
    else
        error('Wrong script to run on this Comp architecture. I cannot detect any windows system')
    end
    fprintf('Mex`s compiled correctly\n')

repo: sr-metric


text:
classdef RWMOP50 < PROBLEM
% <multi> <real> <constrained>
% Power distribution system planning

%------------------------------- Reference --------------------------------
% A. Kumar, G. Wu, M. Ali, Q. Luo, R. Mallipeddi, P. Suganthan, and S. Das,
% A benchmark-suite of real-world constrained multi-objective optimization
% problems and some baseline results, Swarm and Evolutionary Computation,
% 2021, 67: 100961.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Initialization
        function Setting(obj)
            obj.M        = 2;
            obj.D        = 6;
            obj.lower    = [10 10 35 35 125 130];
            obj.upper    = [125 150 210 225 315 325];
            obj.encoding = ones(1,obj.D);
        end
        %% Evaluate multiple solutions
        function Population = Evaluation(obj,varargin)
            x = varargin{1};
            PD = 1200;
            B = [140 17 15 19 26 22;
                  17 60 13 16 15 20;
                  15 13 65 17 24 19;
                  19 16 17 71 30 25;
                  26 15 24 30 69 32;
                  22 20 19 25 32 85] * 10^-6;
            a     = [756.7988 451.3251 1243.5311 1049.9977 1356.6592 1658.5696];
            b     = [38.5390  46.1591 38.3055  40.3965 38.2704 36.3278];
            c     = [0.15247  0.10587  0.03546  0.02803  0.01799 0.02111];
            alpha = [13.8593 13.8593  40.2669  40.2669 42.8955 42.8955];
            beta  = [0.32767  0.32767 -0.54551 -0.54551 -0.51116 -0.51116];
            gamma = [0.00419 0.00419  0.00683  0.00683 0.00461 0.00461];
            PL = zeros(size(x,1),1);
            for i = 1 : size(x,2)
                for j = 1 : size(x,2)
                    PL = PL + x(:,i) .* B(i,j) .* x(j);
                end
            end
            % Objectives
            f = zeros(size(x,1),2);
            for i = 1 : size(x,2)
                f(:,1) = f(:,1) + a(i) + b(i) .* x(:,i) + c(i) .* x(:,i).^2 ;
            end
            for i = 1 : size(x,2)
                f(:,2) = f(:,2) + alpha(i) + beta(i) .* x(:,i) + gamma(i) .* x(:,i).^2;
            end
            % Constraints
            h = sum(x,2) - PD - PL;
            Population = SOLUTION(varargin{1},f,abs(h)-1e-4,varargin{2:end});
            obj.FE     = obj.FE + length(Population);
        end
        %% Generate a point for hypervolume calculation
        function R = GetOptimum(obj,~)
            R = [6.4895619e+04   1.2538936e+03];
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function net = net_init(opts)
% CNN_MNIST_LENET Initialize a CNN similar for MNIST

rng('default');
rng(0) ;

f=1/100 ;
net.layers = {} ;
net.layers{end+1} = struct('type', 'linear', ...
                           'weights', {{f*randn(128,28*28*1, 'single'), zeros(128,1,'single')}}) ;
net.layers{end+1} = struct('type', 'relu') ;
%
net.layers{end+1} = struct('type', 'linear', ...
                           'weights', {{f*randn(128,128, 'single'),  zeros(128,1,'single')}}) ;
net.layers{end+1} = struct('type', 'relu') ;

%
%}
net.layers{end+1} = struct('type', 'linear', ...
                           'weights', {{f*randn(10,128, 'single'), zeros(10,1,'single')}}) ;

net.layers{end+1} = struct('type', 'softmaxloss') ;

for i=1:numel(net.layers)
    if strcmp(net.layers{i}.type,'linear')
        net.layers{1,i}.momentum{1}=zeros(size(net.layers{1,i}.weights{1}));
        net.layers{1,i}.momentum{2}=zeros(size(net.layers{1,i}.weights{2}));
    end
end

%}

author: yechengxi

repo: LightNet

save_path: github-repos/MATLAB/yechengxi-LightNet


text:
function [nb_sources]=source_number_detection_IC(received_signal,criterion,display)

%HELP source_number_detection_IC
%
%Detection of the number of sources with information criteria. The method
%assumes that the noise is spatially white.
%
%Input:     received_signal
%           criterion ('Akaike' or 'MDL')
%           (optional) Display criterion if non empty
%Sorties:   nb_sources
%
%Example:   X=randn(3,500); %3 sources
%           H=randn(10,3)+i*randn(10,3); %10 receivers
%           B=2*(randn(10,500)+i*randn(10,500));
%           Y=H*X+B; %MIMO Mixing Model
%           source_number_detection_IC(Y,'MDL',1)
%
%Reference: [WAX85] Wax, M. and Kailath, T., "Detection of signals by information
%theoretic criteria", IEEE Transactions on Acoustics, Speech and Signal
%Processing, 1985.

[Nb_receivers,Nb_samples]=size(received_signal);
%Computation of the covariance matrix
covariance=received_signal*received_signal'/Nb_samples;
%eigenvalue decomposition
[U,V]=eig(covariance);
%sort the eigenvalue in the decreasing order
eigenvalues=sort(diag(V),'descend');

%Compute criterion
for k=0:Nb_receivers-1
    coef=1/(Nb_receivers-k);
    a=coef*sum(eigenvalues(k+1:Nb_receivers));
    g=prod(eigenvalues(k+1:Nb_receivers)).^(coef);
    akaike_criterion(k+1)=-log(((g/a)^(Nb_samples*(Nb_receivers-k))))+k*(2*Nb_receivers-k);
    MDL_criterion(k+1)=-log(((g/a)^(Nb_samples*(Nb_receivers-k))))+0.5*k*(2*Nb_receivers-k)*log(Nb_samples);
end

switch criterion
    case 'AIC'
        criterion_value=akaike_criterion;
    case 'MDL'
        criterion_value=MDL_criterion;
end
[criterion_value_min,nb_sources]=min(criterion_value);
nb_sources=nb_sources-1; %retrieve 1 because the value k=0 is stocked in the index 1.

if nargin > 2
   plot(0:Nb_receivers-1,criterion_value);
   grid;
   hold on
   plot(nb_sources,criterion_value_min,'ro');
   hold off
   ylabel('Criterion value')
   xlabel('Number of sources')
end

%Created by: Vincent Choqueuse, PhD (vincent.choqueuse@gmail.com)


text:
function a = gfpp_inverse ( n, alpha )

%*****************************************************************************80
%
%% GFPP_INVERSE returns the inverse of the GFPP matrix.
%
%  Example:
%
%    N = 5, ALPHA = 1
%
%    0.5000   -0.2500   -0.1250   -0.0625   -0.0625
%         0    0.5000   -0.2500   -0.1250   -0.1250
%         0         0    0.5000   -0.2500   -0.2500
%         0         0         0    0.5000   -0.5000
%    0.5000    0.2500    0.1250    0.0625    0.0625
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 April 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real ALPHA, determines subdiagonal elements.
%
%    Output, real A(N,N), the inverse matrix.
%
  [ p, l, u ] = gfpp_plu ( n, alpha );

  p_inverse = p';

  l_inverse = tri_l1_inverse ( n, l );

  u_inverse = tri_u_inverse ( n, u );

  a = u_inverse * l_inverse * p_inverse;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function contourdemo

% Demo function for tricontour
%
% Darren Engwirda - 2006

clc, close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                        Driven Cavity
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

answer = lower(input(['This is a demo function for tricontour. \n'                                 ...
                      '\n'                                                                         ...
                      'The following meshes were generated using my mesh generator, "mesh2d.m" \n' ...
                      'and the data comes from my CFD code "Navier2d.m" \n'                        ...
                      '\n'                                                                         ...
                      'Continue?? [y/n] \n'],'s'));

if ~strcmp(answer,'y')
    return
end

load driven_cavity.mat

node  = old_data.node;
cnect = old_data.cnect;
p     = old_data.p;
t     = old_data.t;
U     = old_data.U;
V     = old_data.V;
P     = old_data.P;

figure

subplot(1,2,1), trimesh(t,p(:,1),p(:,2),U), axis square, title('X velocity')
subplot(1,2,2), trimesh(t,p(:,1),p(:,2),V), axis square, title('Y velocity')

figure

set(gcf,'Name','The xy velocity components for a box flow')

subplot(1,2,1), tricontour(p,t,U,15); axis equal, axis off, title('X velocity')
subplot(1,2,2), tricontour(p,t,V,15); axis equal, axis off, title('Y velocity')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                          Tester
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

answer = lower(input(['As well as making pretty pictures, the correct contouring interval is \n'  ...
                      'actually represented. \n'                                                  ...
                      '\n'                                                                        ...
                      'The following shows the contours of distance from a point at [0.5,0.5] \n' ...
                      'The function is correctly showing the contours at a distance of \n'        ...
                      '[0.1,0.2,0.3,0.4] \n'                                                      ...
                      '\n'                                                                        ...
                      'Continue?? [y/n] \n'],'s'));

if ~strcmp(answer,'y')
    return
end

xc = 0.5;
yc = 0.5;
d  = sqrt( (p(:,1)-xc).^2+(p(:,2)-yc).^2 );

figure, trimesh(t,p(:,1),p(:,2),d)
figure, [c,h]=tricontour(p,t,d,[0.1,0.2,0.3,0.4]); clabel(c,h), axis equal, grid on, hold on, plot(0.5,0.5,'bx')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                       Vortex Shedding
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

answer = lower(input(['The domain can be complex. The following example is a domain with holes. \n'             ...
                      '\n'                                                                                ...
                      'The complexity of the domain shouldn''t matter, you just need a valid triangulation. \n' ...
                      '\n'                                                                                ...
                      'Continue?? [y/n] \n'],'s'));

if ~strcmp(answer,'y')
    return
end

close all

load vortex_shedding.mat

node  = old_data.node;
cnect = old_data.cnect;
p     = old_data.p;
t     = old_data.t;
U     = old_data.U;
V     = old_data.V;
P     = old_data.P;
W     = old_data.W;

figure

trimesh(t,p(:,1),p(:,2),U), axis square, title('X velocity')

figure

set(gcf,'Name','The x velocity for the flow over a cylinder')

% Velocity contours
tricontour(p,t,U,50), axis equal, axis off, title('X velocity')

% Walls
patch('faces',cnect,'vertices',node,'facecolor','none','edgecolor','k')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                        Peaks
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

answer = lower(input(['A standard contouring benchmark showing the operation of clabel. \n' ...
                      '\n'                                                                  ...                                                                                ...
                      'Continue?? [y/n] \n'],'s'));

if ~strcmp(answer,'y')
    return
end

close all

[xx,yy] = meshgrid(linspace(-3,3,64),linspace(-2.5,2.5,64));
zz      = peaks(xx,yy);
v       = -3:5;
figure(1)
[c,h] = contour(xx,yy,zz,v);
clabel(c,h)
title('Contour')

% Triangulate
p = [xx(:),yy(:)];
t = delaunayn(p);

figure(2)
[c,h] = tricontour(p,t,zz(:),v);
clabel(c,h)
title('Tricontour')


text:
function test_bug2647
% DEPENDENCY ft_freqstatistics

% WALLTIME 00:10:00
% MEM 2gb

load(dccnpath('/home/common/matlab/fieldtrip/data/test/bug2647.mat'));
stat = ft_freqstatistics(cfg,data,BL);
%assert(all(~isfinite(stat.stat(:))));
assert(~all(~isfinite(stat.stat(:)))); % this should now work after making ft_statfun_actvsblT more robust for NaNs

cfg.design = cfg.design(:,[1:4 6:9]);
data.powspctrm = data.powspctrm(2:end,:,:,:);
BL.powspctrm   = BL.powspctrm(2:end,:,:,:);
stat = ft_freqstatistics(cfg,data,BL);
assert(~all(~isfinite(stat.stat(:))));

% conclusion: the first rpt was all NaN, causing stat.stat to be all NaN,
% no meaningful inference possible

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function host = gethostname()

% HOSTNAME returns the hostname of this computer
%
% Use as
%   str = hostname;
%
% See also GETUSERNAME, GETADDRESS

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

% this is to speed up subsequent calls
persistent previous_argout
if ~isempty(previous_argout)
  host = previous_argout;
  return
end

if (ispc())
  host = getenv('ComputerName');
elseif (isunix())
  host = getenv('HOSTNAME');
  % the HOSTNAME variable is not guaranteed to be set
  if isempty(host)
    [status, host] = system('hostname -s');
    % remove the ^n from the end of the line
    host = deblank(host);
  end
end

host = strtok(host, '.'); % dots in filenames are not allowed by matlab

if (isempty(host))
  host = 'unknownhost';
end

% remember for subsequent calls
previous_argout = host;


text:
function display_bbox_detections( img, bbox_detections, score_thresh, category_names )

all_dets = [];
num_categories = length(bbox_detections);
for i = 1:num_categories
    bbox_detections_this_category  = bbox_detections{i};
    detection_scores_this_category = bbox_detections_this_category(:,5);
    is_above_the_thresh = detection_scores_this_category >= score_thresh(i);
    bbox_detections_this_category = bbox_detections_this_category(is_above_the_thresh,:);
    bbox_detections_this_category = [i * ones(size(bbox_detections_this_category, 1), 1), bbox_detections_this_category];
    all_dets = cat(1, all_dets, bbox_detections_this_category);
end

fprintf('Visualize the bounding box detections:\n')
[~, ord] = sort(all_dets(:,end), 'descend');
for i = 1:length(ord)
  score_this = all_dets(ord(i), end);
  category_name_this = category_names{all_dets(ord(i), 1)};
  showboxes(img, all_dets(ord(i), 2:5));
  title(sprintf('det #%d: %s score = %.3f', i, category_name_this, score_this));
  fprintf('det #%d: %s score = %.3f. press any key to continue\n', ...
      i, category_name_this, score_this);
  drawnow;
  pause;
end

fprintf('No more detections\n');

end

author: gidariss

repo: LocNet

save_path: github-repos/MATLAB/gidariss-LocNet


text:
function intVal=ellipIntInc1Kind(phi,m)
%%ELLIPINTINC1KIND Evaluate Legendre's incomplete elliptical integral of
%                  the first kind. This is the integral from 0 to phi of
%                  (1-m*sin(theta)^2)^(-1/2) dtheta. Alternatively, this
%                  can be written as the integral from 0 to x of
%                  ((1-t^2)*(1-m*t^2))^(-1/2) dt where x=sin(phi). This
%                  integral is sometimes referred to as F.
%
%INPUTS: phi A scalar or matrix of the real upper bounds of integration.
%            These can be positive, negative or zero.
%          m A value between 0 and 1. The integral diverges when m=1 and
%            phi=pi/2.
%
%OUTPUTS: intVal The values of the incomplete elliptic integral of the
%                first kind evaluated at each entry in phi.
%
%The expression used to implement the function comes from [1], whereby
%negative values are handled by switching the sign of the result. The above
%algorithm also only considers values with magnitudes up to pi/2. To get
%the correct result, however, for larger values, the value up to pi/2
%(times the number of multiples of the value up to pi/2) is added to a
%fractional part with magnitude less than pi/2. Note however, that due to
%the odd symmetry of sin(x)^2, if the integer part is odd, then the
%computation subtracts the fractional part from pi/2 and the value added is
%the value of the complete integral MINUS the modified fractional part.
%
%REFERENCES:
%[1] B. C. Carlson, "Numerical computation of real or complex elliptic
%    integrals," Numerical Algorithms, vol. 10, no. 1, pp. 13-26, 1995.
%
%September 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numItems=length(phi(:));
intVal=zeros(size(phi));

for curItem=1:numItems
    phiCur=phi(curItem);

    phiSign=sign(phiCur);
    phiCur=abs(phiCur);

    %Separate out the multiples of pi/2 from the rest of the values.
    intPart=fix(phiCur/(pi/2));
    if(intPart>0)
        completeVal=symIntFirstKind(0,1-m,1);
        phiCur=phiCur-intPart*(pi/2);
    else
        completeVal=0;
    end

    %The test deals with precision limitations that would prevent the other
    %integral from converging.
    if(phiCur<eps)
        intVal(curItem)=phiSign*intPart*completeVal;
    else
        %Deal with the odd symmetry of sin(theta)^2
        if(mod(intPart,2)==0)
            c=csc(phiCur)^2;
            delta=symIntFirstKind(c-1,c-m,c);

            intVal(curItem)=phiSign*(intPart*completeVal+delta);
        else
            phiCur=pi/2-phiCur;

            c=csc(phiCur)^2;
            delta=symIntFirstKind(c-1,c-m,c);

            intVal(curItem)=phiSign*(intPart*completeVal+(completeVal-delta));
        end
    end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function r8_big_test ( )

%*****************************************************************************80
%
%% R8_BIG_TEST tests R8_BIG.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 November 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8_BIG_TEST\n' );
  fprintf ( 1, '  R8_BIG returns a "big" R8 value;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  R8_BIG = %f\n', r8_big ( ) );

  return
end


text:
% SIGNED_DISTANCE Compute signed distance from points P to a mesh (V,F)
%
% [S,I,C,N] = signed_distance(P,V,F,'ParameterName',parameter_value,...)
%
% Inputs:
%   P  #P by 3 list of query point positions
%   V  #V by 3 list of vertex positions
%   F  #F by 3 list of triangle indices
%   Optional:
%     'SignedDistanceType' followed by
%       'winding_number' use winding number (continuous sign value for
%         non-watertight)
%       {'pseudonormal'}  use pseudo-normal, binary scale (but not robust for
%         non-watertight meshes.
% Outputs:
%   S  #P list of smallest signed distances
%   I  #P list of facet indices corresponding to smallest distances
%   C  #P by 3 list of closest points
%   N  #P by 3 list of closest normals (only set if
%
% Example:
%

author: alecjacobson

repo: gptoolbox


text:
function mgdata = mg_ns_diff(Ap,domain,flowsol,qmethod,outbnd)
%mg_ns_diff   GMG for diffusion problem (Navier-Stokes)
%   mgdata = mg_ns_diff(Ap,domain,flowsol,qmethod,outbnd)
%   input
%          Ap           top level discrete diffusion operator
%          domain       domain information, 1 for cavity, 3 for step
%          flowsol      current velocity solution iterate
%          qmethod      discretization method, 0 for Q1-Q1, 2 for Q2-Q1
%          outbnd       location of outflow boundary (for step)
%   output
%          mgdata       structure containing GMG data at all levels
%             matrix    discrete diffusion operators
%             prolong   grid transfer operators
%             smoother  structure containing smoothing operators
%             nc        grid parameter identifying finest level
%
%   IFISS function: HCE; 18 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage

if domain==1,      % Cavity
   np = sqrt(length(Ap));
   nc = log2(np-1);
   grid_type = 1;

% Smoothing strategy
   sweeps=1; smooth=2; stype=2;  % Means one (sweeps)-directional line (stype) GS (smooth)

% Top level
% Note (19 April 2004):  Fine level matrix created here is not pinned down
   [x,y,xy,mv,mp,bound] = mg_cavity_domain(nc+1,grid_type);
   [x,y,xy,xyp] = mg_diff_q2q1grid(x,y,xy,mv,mp,bound);
   mgdata(nc).matrix = mg_apsetup_q1(xy,xyp,mv,mp,domain);
   if qmethod==2,
      [x,y,xy,mv,mp,bound] = mg_cavity_domain(nc,grid_type);
      [x,y,xy,xyp] = mg_diff_q2q1grid(x,y,xy,mv,mp,bound);
      mgdata(nc).prolong = mg_prolong(2^nc,2^nc,x,y);
   else % qmethod=0
      mgdata(nc).prolong=mg_prolong(2^nc,2^nc,x,y);
   end
   mgdata(nc).smoother = mg_ns_smooth(mgdata(nc).matrix,sweeps,smooth,stype);
   mgdata(nc).nc = nc;  % Not used, saved to be consistent with domain=3

   % Loop over remaining levels
   for level = nc-1:-1:2;
      if qmethod==2,
         %%% construct matrix
         mgdata(level).matrix = mg_apsetup_q1(xy,xyp,mv,mp,domain);
         %%% construct prolongation operator
         [x,y,xy,mv,mp,bound] = mg_cavity_domain(level,grid_type);
         [x,y,xy,xyp] = mg_diff_q2q1grid(x,y,xy,mv,mp,bound);
         mgdata(level).prolong = mg_prolong(2^level,2^level,x,y);
      else % qmethod=0
         [x,y,xy,mv,mp,bound] = mg_cavity_domain(level+1,grid_type);
         [x,y,xy,xyp] = mg_diff_q2q1grid(x,y,xy,mv,mp,bound);
         mgdata(level).matrix = mg_apsetup_q1(xy,xyp,mv,mp,domain);
         mgdata(level).prolong=mg_prolong(2^level,2^level,x,y);
      end
      mgdata(level).smoother = mg_ns_smooth(mgdata(level).matrix,sweeps,smooth,stype);
      mgdata(level).nc = nc;
   end

elseif domain==3,     % Step
   nu = length(flowsol)/2;
   xi = (-(outbnd+3)+sqrt((outbnd+3)^2+4*(2*outbnd+1)*(nu-1))) ...
          / (2*(2*outbnd+1));
   nc = log2(2*xi);

% Smoothing strategy
   sweeps=1; smooth=3; stype=1; % Means one (sweeps)-directional point (stype) ILU (smooth)

% top level
   [x,y,xy,mv,mp,bound] = mg_step_domain(nc,outbnd);
   mgdata(nc).matrix = Ap;
   if qmethod==2,
      x=x(1:2:end)'; y=y(1:2:end)'; mgdata(nc).prolong=mg_ns_prolong_step(nc-1,x,y,outbnd);
   else % qmethod=0
      x=x'; y=y';
      mgdata(nc).prolong=mg_ns_prolong_step(nc,x,y,outbnd);
   end
   mgdata(nc).smoother = mg_ns_smooth(mgdata(nc).matrix,sweeps,smooth,stype);
   mgdata(nc).nc = nc;

%%% loop over remaining levels
   for level = nc-1:-1:2;
      if qmethod==2,
         [x,y,xy,mv,mp,bound] = mg_step_domain(level,outbnd);
         [x,y,xy,xyp] = mg_diff_q2q1grid(x,y,xy,mv,mp,bound);
         mgdata(level).matrix = mg_apsetup_q1(xy,xyp,mv,mp,domain);
         if level>2,
            x=x(1:2:end)'; y=y(1:2:end)';
            mgdata(level).prolong=mg_ns_prolong_step(level-1,x,y,outbnd);
         end
      else % qmethod=0
         [x,y,xy,xyp] = mg_diff_q2q1grid(x,y,xy,mv,mp,bound);
         mgdata(level).matrix = mg_apsetup_q1(xy,xyp,mv,mp,domain);
         [x,y,xy,mv,mp,bound] = mg_step_domain(level,outbnd);
         x=x'; y=y';
         mgdata(level).prolong = mg_ns_prolong_step(level,x,y,outbnd);
      end
      mgdata(level).smoother = mg_ns_smooth(mgdata(level).matrix,sweeps,smooth,stype);
      mgdata(level).nc = nc;
   end
else
   error('Error, mg_ns_diff_data is defined only for cavity and step\n');
end


text:
%-------------------------These are the input parameters----------------%
Input_Depth     = 'Ballet40.yuv';       %This is the input depth video
Output_Depth    = 'Output_Ballet.yuv';  %Give a name to the output file
Color           = 'Ballet.yuv';         %Corresponding color video
Width           = 1024;                 %Width in pixels
Height          = 768;                  %Height in pixels
Frame_Num       = 1;                    %Number of Frames to process
%------------------------End of the Input Parameters--------------------%

%-------------------------Key Performance Parameters--------------------%
w               = 15;                   %Half width of the kernel (odd)
sigma_range     = 0.025                 %Standard Deviation of the range
                                        %filter, change this suitably.
%------------------------End of key perf. parameters--------------------%

sigma = [3*w sigma_range]; % bilateral filter standard deviations

% Apply Joint Bilateral Filter
filename_in = Input_Depth;
fid_in = fopen(filename_in, 'rb');
[Yd, Ud, Vd] = yuv_import(filename_in,[Width Height],Frame_Num);
fclose(fid_in);

filename_in = Color;
fid_in = fopen(filename_in, 'rb');
[Yc, Uc, Vc] = yuv_import(filename_in,[Width Height],Frame_Num);
fclose(fid_in);

fid = fopen(Output_Depth, 'w');

for s = 1:Frame_Num,
    Ad = Yd{s};
    Ac = Yc{s};
    Ax = jbfilter2(double(Ad)/255,double(Ac)/255,w,sigma);
    image = round(Ax*255);
    FramesOut.Luma = image;
    FramesOut.Chroma1 = Ud{s};
    FramesOut.Chroma2 = Vd{s};
    figure, imshow(uint8(FramesOut.Luma));
    fwrite(fid, FramesOut.Luma', 'uint8');
    fwrite(fid, FramesOut.Chroma1', 'uint8');
    fwrite(fid, FramesOut.Chroma2', 'uint8');
end

fclose(fid);


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% Tutorial for FAIR:  Landmark Based Registration, quadratic transformation
%
% - load data (see setup2DhandData)
% - setup  viewer (viewImage2D), interpolator (splineInter),
% - setup landmarks (LM)
% - run quadratic
%==============================================================================

clear, close all, help(mfilename)

%% setup hand data
setup2DhandData

if FAIRinput('','set new landmarks ? ',0),
  [LM,fig] = getLandmarks(dataT,dataR,omega,m);
  close(fig);
end;

omegaT = omega(1,:);
omegaR = omega(end,:);
xT = getCellCenteredGrid(omegaT,m);
xR = getCellCenteredGrid(omegaR,m);
Tc = imgModel(dataT,omegaT,xT);
Rc = imgModel(dataR,omegaR,xR);

%% visualize data
FAIRfigure(1,'figname',mfilename); clf;
subplot(1,3,1); viewImage(Tc,omegaT,m); hold on;
ph = plotLM(LM(:,1:2),'numbering','on','color','r');
set(ph,'linewidth',2,'markersize',20);
title(sprintf('%s','T&LM'),'fontsize',20);

subplot(1,3,2); viewImage(Rc,omegaR,m); hold on;
ph = plotLM(LM(:,3:4),'numbering','on','color','g','marker','+');
set(ph,'linewidth',2,'markersize',20);
title(sprintf('%s','R&LM'),'fontsize',20);

%% compute landmark based registration
[yc,LM] = LMreg('quadratic',LM(:,1:4),xR);
TLM = imgModel(dataT,omegaT,yc);

subplot(1,3,3); cla; viewImage(TLM,omegaR,m); hold on;
ph = plotLM(LM(:,3:4),'numbering','off','color','g','marker','+');
qh = plotLM(LM(:,7:8),'numbering','off','color','m','marker','x');
rh = plot(LM(:,[3,7])',LM(:,[4,8])','m-','linewidth',3);
set([ph;qh;rh],'linewidth',2,'markersize',20);
title(sprintf('%s','T(y^{quadratic})&LM'),'fontsize',20);
%==============================================================================


text:
function combo_test27 ( )

%*****************************************************************************80
%
%% COMBO_TEST27 tests PERM_INV and PERM_MUL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2011
%
%  Author:
%
%    John Burkardt
%
  n = 4;

  fprintf ( 1, ' \n' );
  fprintf ( 1, 'COMBO_TEST27\n' );
  fprintf ( 1, '  Permutations of the integers:\n' );
  fprintf ( 1, ' \n' );
  fprintf ( 1, '  PERM_INV computes an inverse permutation,\n' );
  fprintf ( 1, '  PERM_MUL multiplies two permutations.\n' );
%
%  Enumerate.
%
  nperm = perm_enum ( n );
%
%  Unrank.
%
  rank = floor ( nperm / 2 );

  p = perm_lex_unrank ( rank, n );

  perm_print ( n, p, '  The permutation P is ' );
%
%  Invert.
%
  q = perm_inv ( n, p );

  perm_print ( n, q, '  The inverse permutation Q is ' );
%
%  Multiply.
%
  r = perm_mul ( n, p, q );

  perm_print ( n, r, '  The product R = P * Q is ' );

  return
end


text:
classdef PTKMetrics < dynamicprops
    % PTKMetrics. A structure for holding data analysis results
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties (SetAccess = private)
        MetricNameMap
    end

    methods
        function obj = PTKMetrics
            obj.MetricNameMap = containers.Map;
        end

        function AddMetric(obj, property_name, value, user_visible_name)
            valid_property_name = CoreTextUtilities.CreateValidFieldName(property_name);
            obj.addprop(valid_property_name);
            obj.(valid_property_name) = value;
            obj.MetricNameMap(valid_property_name) = user_visible_name;
        end

        function property_list = GetListOfMetrics(obj)
            property_list = properties(obj);
            property_list = setdiff(property_list, 'MetricNameMap');
        end

        function Merge(obj, metrics, reporting, prefix, user_visible_prefix)
            if nargin < 4
                prefix = [];
                user_visible_prefix = [];
            end
            if ~isa(metrics, 'PTKMetrics')
                reporting.Error('PTKMetrics:NotAMetricsClass', 'The argument passed to Merge() was not of type PTKMetrics');
            end

            metric_list = metrics.GetListOfMetrics;
            for i = 1 : length(metric_list)
                metric = metric_list{i};
                obj.AddMetric([prefix, metric], metrics.(metric), [user_visible_prefix, metrics.MetricNameMap(metric)]);
            end
        end
    end

    methods (Static)
        function results = MergeResults(results_1, results_2)
            if isempty(results_2)
                results = results_1;
            else
                metric_names = unique([fieldnames(results_1), fieldnames(results_2)]);
                results = [];
                for metric_cell = metric_names'
                    metric_name = metric_cell{1};
                    if isfield(results_1, metric_name)
                        results.(metric_name) = results_1.(metric_name);
                    end
                    if isfield(results, metric_name)
                        if isfield(results_2, metric_name)
                            results.(metric_name).Merge(results_2.(metric_name))
                        else
                            results.(metric_name) = results_2.(metric_name);
                        end
                    end
                end
            end
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function a = imag(a)
%IMAG         imaginary part of gradients
%
%   c = imag(a)
%
%

% written  10/16/98     S.M. Rump
% modified 12/18/02     S.M. Rump  improved performance
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  a.x = imag(a.x);
  a.dx = imag(a.dx);

repo: OpenMAS


text:
%DEMO_SVI_REGRESSION  A toy data regression example for sparse SVI GP
%
%  Description
%    Demonstration of sochastic variational inference GP model regression.
%    The problem is similar as in the Hensman et. al (2013). The dataset is
%    made synthetically in two dimensions using sinusoidal functions.
%
%  See also
%    DEMO_SVI_CLASSIFIC
%
%  References:
%    Hensman, J., Fusi, N. and Lawrence, N. D. (2013). Gaussian processes
%    for big data. arXiv preprint arXiv:1309.6835.

% Copyright (c) 2014 Tuomas Sivula

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

% Generate toy data
n = 2000;
s2 = 0.01;
% Training samples
x = rand(n,2)*2*pi - pi;
y = sin(x(:,1)).*sin(x(:,2)) + sqrt(s2).*randn(n,1);
% Test samples
nt = 100;
xt = rand(nt,2)*2*pi - pi;
yt = sin(xt(:,1)).*sin(xt(:,2));

% Prediction grid
[X1,X2] = meshgrid(-4:0.2:4, -4:0.2:4);
grid = [X1(:), X2(:)];

% Build gp structure
gp = gp_set('lik',lik_gaussian, 'cf', gpcf_sexp, 'latent_method', 'SVI');
% Optimise
nu = 50;    % The number of inducing inputs
maxi = 150; % The maximum number of iteration rounds
[gp, diagnosis] = svigp(gp,x,y,'xt',xt,'yt',yt,'nu',nu,'maxiter',maxi);
% Predict
Eft = gpsvi_pred(gp,x,y,grid);

% ------- Plot --------
% Prediction
figure()
contour(X1, X2, reshape(Eft, size(X1)), 16)
hold on
scatter(x(:,1), x(:,2), 20, y, 'filled')
scatter(gp.X_u(:,1), gp.X_u(:,2), 30)
legend('Eft', 'data', 'Z')
% Convergence analysis
figure()
subplot(3,1,1)
plot(mean(diagnosis.e,2))
title('energy')
subplot(3,1,2)
plot(diagnosis.mlpd)
title('mean log predictive density')
subplot(3,1,3)
plot(diagnosis.rmse)
title('root mean square error')
xlabel('iteration')

author: gpstuff-dev

repo: gpstuff

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