text:
function pesq_mos= pesq_psychoacoustic_model (ref_data, ref_Nsamples, deg_data, ...
    deg_Nsamples )

global CALIBRATE Nfmax Nb Sl Sp
global nr_of_hz_bands_per_bark_band centre_of_band_bark
global width_of_band_hz centre_of_band_hz width_of_band_bark
global pow_dens_correction_factor abs_thresh_power
global Downsample SEARCHBUFFER DATAPADDING_MSECS Fs Nutterances
global Utt_Start Utt_End Utt_Delay NUMBER_OF_PSQM_FRAMES_PER_SYLLABE
global Fs Plot_Frame

% Plot_Frame= 75; % this is the frame whose spectrum will be plotted

FALSE= 0;
TRUE= 1;
NUMBER_OF_PSQM_FRAMES_PER_SYLLABE= 20;

maxNsamples = max (ref_Nsamples, deg_Nsamples);
Nf = Downsample * 8;
MAX_NUMBER_OF_BAD_INTERVALS = 1000;

start_frame_of_bad_interval= zeros( 1, MAX_NUMBER_OF_BAD_INTERVALS);
stop_frame_of_bad_interval= zeros( 1, MAX_NUMBER_OF_BAD_INTERVALS);
start_sample_of_bad_interval= zeros( 1, MAX_NUMBER_OF_BAD_INTERVALS);
stop_sample_of_bad_interval= zeros( 1, MAX_NUMBER_OF_BAD_INTERVALS);
number_of_samples_in_bad_interval= zeros( 1, MAX_NUMBER_OF_BAD_INTERVALS);
delay_in_samples_in_bad_interval= zeros( 1, MAX_NUMBER_OF_BAD_INTERVALS);
number_of_bad_intervals= 0;
there_is_a_bad_frame= FALSE;

Whanning= hann( Nf, 'periodic');
Whanning= Whanning';

D_POW_F = 2;
D_POW_S = 6;
D_POW_T = 2;
A_POW_F = 1;
A_POW_S = 6;
A_POW_T = 2;
D_WEIGHT= 0.1;
A_WEIGHT= 0.0309;

CRITERIUM_FOR_SILENCE_OF_5_SAMPLES = 500;
samples_to_skip_at_start = 0;
sum_of_5_samples= 0;
while ((sum_of_5_samples< CRITERIUM_FOR_SILENCE_OF_5_SAMPLES) ...
        && (samples_to_skip_at_start < maxNsamples / 2))
    sum_of_5_samples= sum( abs( ref_data( samples_to_skip_at_start...
        + SEARCHBUFFER * Downsample + 1: samples_to_skip_at_start...
        + SEARCHBUFFER * Downsample + 5)));

    if (sum_of_5_samples< CRITERIUM_FOR_SILENCE_OF_5_SAMPLES)
        samples_to_skip_at_start = samples_to_skip_at_start+ 1;
    end
end
% fprintf( 'samples_to_skip_at_start is %d\n', samples_to_skip_at_start);

samples_to_skip_at_end = 0;
sum_of_5_samples= 0;
while ((sum_of_5_samples< CRITERIUM_FOR_SILENCE_OF_5_SAMPLES) ...
        && (samples_to_skip_at_end < maxNsamples / 2))
    sum_of_5_samples= sum( abs( ref_data( maxNsamples - ...
        SEARCHBUFFER* Downsample + DATAPADDING_MSECS* (Fs/ 1000) ...
        - samples_to_skip_at_end - 4: maxNsamples - ...
        SEARCHBUFFER* Downsample + DATAPADDING_MSECS* (Fs/ 1000) ...
        - samples_to_skip_at_end)));
    if (sum_of_5_samples< CRITERIUM_FOR_SILENCE_OF_5_SAMPLES)
        samples_to_skip_at_end = samples_to_skip_at_end+ 1;
    end
end
% fprintf( 'samples_to_skip_at_end is %d\n', samples_to_skip_at_end);

start_frame = floor( samples_to_skip_at_start/ (Nf/ 2));
stop_frame = floor( (maxNsamples- 2* SEARCHBUFFER* Downsample ...
    + DATAPADDING_MSECS* (Fs/ 1000)- samples_to_skip_at_end) ...
    / (Nf/ 2))- 1;
% number of frames in speech data plus DATAPADDING_MSECS
% fprintf( 'start/end frame is %d/%d\n', start_frame, stop_frame);

D_disturbance= zeros( stop_frame+ 1, Nb);
DA_disturbance= zeros( stop_frame+ 1, Nb);

power_ref = pow_of (ref_data, SEARCHBUFFER* Downsample, ...
    maxNsamples- SEARCHBUFFER* Downsample+ DATAPADDING_MSECS* (Fs/ 1000),...
    maxNsamples- 2* SEARCHBUFFER* Downsample+ DATAPADDING_MSECS* (Fs/ 1000));
power_deg = pow_of (deg_data, SEARCHBUFFER * Downsample, ...
    maxNsamples- SEARCHBUFFER* Downsample+ DATAPADDING_MSECS* (Fs/ 1000),...
    maxNsamples- 2* SEARCHBUFFER* Downsample+ DATAPADDING_MSECS* (Fs/ 1000));
% fprintf( 'ref/deg power is %f/%f\n', power_ref, power_deg);

hz_spectrum_ref             = zeros( 1, Nf/ 2);
hz_spectrum_deg             = zeros( 1, Nf/ 2);
frame_is_bad                = zeros( 1, stop_frame + 1);
smeared_frame_is_bad        = zeros( 1, stop_frame + 1);
silent                      = zeros( 1, stop_frame + 1);

pitch_pow_dens_ref          = zeros( stop_frame + 1, Nb);
pitch_pow_dens_deg          = zeros( stop_frame + 1, Nb);

frame_was_skipped           = zeros( 1, stop_frame + 1);
frame_disturbance           = zeros( 1, stop_frame + 1);
frame_disturbance_asym_add  = zeros( 1, stop_frame + 1);

avg_pitch_pow_dens_ref      = zeros( 1, Nb);
avg_pitch_pow_dens_deg      = zeros( 1, Nb);
loudness_dens_ref           = zeros( 1, Nb);
loudness_dens_deg           = zeros( 1, Nb);
deadzone                    = zeros( 1, Nb);
disturbance_dens            = zeros( 1, Nb);
disturbance_dens_asym_add   = zeros( 1, Nb);

time_weight                 = zeros( 1, stop_frame + 1);
total_power_ref             = zeros( 1, stop_frame + 1);

% fid= fopen( 'tmp_mat.txt', 'wt');

for frame = 0: stop_frame
    start_sample_ref = 1+ SEARCHBUFFER * Downsample + frame* (Nf/ 2);
    hz_spectrum_ref= short_term_fft (Nf, ref_data, Whanning, ...
        start_sample_ref);

    utt = Nutterances;
    while ((utt >= 1) && ((Utt_Start(utt)- 1)* Downsample+ 1 ...
            > start_sample_ref))
        utt= utt - 1;
    end

    if (utt >= 1)
        delay = Utt_Delay(utt);
    else
        delay = Utt_Delay(1);
    end

    start_sample_deg = start_sample_ref + delay;

    if ((start_sample_deg > 0) && (start_sample_deg + Nf- 1 < ...
            maxNsamples+ DATAPADDING_MSECS* (Fs/ 1000)))
        hz_spectrum_deg= short_term_fft (Nf, deg_data, Whanning, ...
            start_sample_deg);
    else
        hz_spectrum_deg( 1: Nf/ 2)= 0;
    end

    pitch_pow_dens_ref( frame+ 1, :)= freq_warping (...
        hz_spectrum_ref, Nb, frame);
    %peak = maximum_of (pitch_pow_dens_ref, 0, Nb);
    pitch_pow_dens_deg( frame+ 1, :)= freq_warping (...
        hz_spectrum_deg, Nb, frame);

    total_audible_pow_ref = total_audible (frame, pitch_pow_dens_ref, 1E2);
    total_audible_pow_deg = total_audible (frame, pitch_pow_dens_deg, 1E2);
    silent(frame+ 1) = (total_audible_pow_ref < 1E7);

end
% fclose( fid);

avg_pitch_pow_dens_ref= time_avg_audible_of (stop_frame + 1, ...
    silent, pitch_pow_dens_ref, floor((maxNsamples- 2* SEARCHBUFFER* ...
    Downsample+ DATAPADDING_MSECS* (Fs/ 1000))/ (Nf / 2))- 1);
avg_pitch_pow_dens_deg= time_avg_audible_of (stop_frame + 1, ...
    silent, pitch_pow_dens_deg, floor((maxNsamples- 2* SEARCHBUFFER* ...
    Downsample+ DATAPADDING_MSECS* (Fs/ 1000))/ (Nf/ 2))- 1);

% fid= fopen( 'tmp_mat.txt', 'wt');
% fprintf( fid, '%f\n', avg_pitch_pow_dens_deg);
% fclose( fid);

if (CALIBRATE== 0)
    pitch_pow_dens_ref= freq_resp_compensation (stop_frame + 1, ...
        pitch_pow_dens_ref, avg_pitch_pow_dens_ref, ...
        avg_pitch_pow_dens_deg, 1000);
    if (Plot_Frame>= 0) % plot pitch_pow_dens_ref
        figure;
        subplot( 1, 2, 1);
        plot( centre_of_band_hz, 10* log10( eps+ ...
            pitch_pow_dens_ref( Plot_Frame+ 1, :)));
        axis( [0 Fs/2 0 95]); %xlabel( 'Hz'); ylabel( 'Db');
        title( 'reference signal bark spectrum with frequency compensation');
        subplot( 1, 2, 2);
        plot( centre_of_band_hz, 10* log10( eps+ ...
            pitch_pow_dens_deg( Plot_Frame+ 1, :)));
        axis( [0 Fs/2 0 95]); %xlabel( 'Hz'); ylabel( 'Db');
        title( 'degraded signal bark spectrum');
    end

end
% tmp1= pitch_pow_dens_ref';

MAX_SCALE = 5.0;
MIN_SCALE = 3e-4;
oldScale = 1;
THRESHOLD_BAD_FRAMES = 30;
for frame = 0: stop_frame

    total_audible_pow_ref = total_audible (frame, pitch_pow_dens_ref, 1);
    total_audible_pow_deg = total_audible (frame, pitch_pow_dens_deg, 1);
    total_power_ref (1+ frame) = total_audible_pow_ref;

    scale = (total_audible_pow_ref + 5e3)/ (total_audible_pow_deg + 5e3);
    if (frame > 0)
        scale = 0.2 * oldScale + 0.8 * scale;
    end
    oldScale = scale;

    if (scale > MAX_SCALE)
        scale = MAX_SCALE;
    elseif (scale < MIN_SCALE)
        scale = MIN_SCALE;
    end

    pitch_pow_dens_deg( 1+ frame, :) = ...
        pitch_pow_dens_deg( 1+ frame, :) * scale;

    if (frame== Plot_Frame)
        figure;
        subplot( 1, 2, 1);
        plot( centre_of_band_hz, 10* log10( eps+ ...
            pitch_pow_dens_ref( Plot_Frame+ 1, :)));
        axis( [0 Fs/2 0 95]); %xlabel( 'Hz'); ylabel( 'Db');
        subplot( 1, 2, 2);
        plot( centre_of_band_hz, 10* log10( eps+ ...
            pitch_pow_dens_deg( Plot_Frame+ 1, :)));
        axis( [0 Fs/2 0 95]); %xlabel( 'Hz'); ylabel( 'Db');
    end

    loudness_dens_ref = intensity_warping_of (frame, pitch_pow_dens_ref);
    loudness_dens_deg = intensity_warping_of (frame, pitch_pow_dens_deg);
    disturbance_dens = loudness_dens_deg - loudness_dens_ref;

    if (frame== Plot_Frame)
        figure;
        subplot( 1, 2, 1);
        plot( centre_of_band_hz, 10* log10( eps+ ...
            loudness_dens_ref));
        axis( [0 Fs/2 0 15]); %xlabel( 'Hz'); ylabel( 'Db');
        title( 'reference signal loudness density');
        subplot( 1, 2, 2);
        plot( centre_of_band_hz, 10* log10( eps+ ...
            loudness_dens_deg));
        axis( [0 Fs/2 0 15]); %xlabel( 'Hz'); ylabel( 'Db');
        title( 'degraded signal loudness density');
    end

    for band =1: Nb
        deadzone (band) = 0.25* min (loudness_dens_deg (band), ...
            loudness_dens_ref (band));
    end

    for band = 1: Nb
        d = disturbance_dens (band);
        m = deadzone (band);

        if (d > m)
            disturbance_dens (band) = disturbance_dens (band)- m;
%             disturbance_dens (band) = d- m;
        else
            if (d < -m)
                disturbance_dens (band) = disturbance_dens (band)+ m;
%                 disturbance_dens (band) = d+ m;
            else
                disturbance_dens (band) = 0;
            end
        end
    end

    if (frame== Plot_Frame)
        figure;
        subplot( 1, 2, 1);
        plot( centre_of_band_hz, disturbance_dens);
        axis( [0 Fs/2 -1 50]); %xlabel( 'Hz'); ylabel( 'Db');
        title( 'disturbance');
    end
    D_disturbance( frame+ 1, :)= disturbance_dens;

    frame_disturbance (1+ frame) = pseudo_Lp (disturbance_dens, D_POW_F);
    if (frame_disturbance (1+ frame) > THRESHOLD_BAD_FRAMES)
        there_is_a_bad_frame = TRUE;
    end

    disturbance_dens= multiply_with_asymmetry_factor (...
        disturbance_dens, frame, pitch_pow_dens_ref, pitch_pow_dens_deg);

    if (frame== Plot_Frame)
        subplot( 1, 2, 2);
        plot( centre_of_band_hz, disturbance_dens);
        axis( [0 Fs/2 -1 50]); %xlabel( 'Hz'); ylabel( 'Db');
        title( 'disturbance after asymmetry processing');
    end
    DA_disturbance( frame+ 1, :)= disturbance_dens;

    frame_disturbance_asym_add (1+ frame) = ...
        pseudo_Lp (disturbance_dens, A_POW_F);
end
% fid= fopen( 'tmp_mat.txt', 'wt');
% fprintf( fid, '%f\n', frame_disturbance);
% fclose( fid);

frame_was_skipped (1: 1+ stop_frame) = FALSE;

for utt = 2: Nutterances
    frame1 = floor (((Utt_Start(utt)- 1- SEARCHBUFFER )* Downsample+ 1+ ...
        Utt_Delay(utt))/ (Nf/ 2));
    j = floor( floor(((Utt_End(utt-1)- 1- SEARCHBUFFER)* Downsample+ 1+ ...
        Utt_Delay(utt-1)))/(Nf/ 2));
    delay_jump = Utt_Delay(utt) - Utt_Delay(utt-1);
    if (frame1 > j)
        frame1 = j;
    elseif (frame1 < 0)
        frame1 = 0;
    end
%     fprintf( 'frame1, j, delay_jump is %d, %d, %d\n', frame1, ...
%         j, delay_jump);

    if (delay_jump < -(Nf/ 2))
        frame2 = floor (((Utt_Start(utt)- 1- SEARCHBUFFER)* Downsample+ 1 ...
            + max (0, abs (delay_jump)))/ (Nf/ 2)) + 1;

        for frame = frame1: frame2
            if (frame < stop_frame)
                frame_was_skipped (1+ frame) = TRUE;
                frame_disturbance (1+ frame) = 0;
                frame_disturbance_asym_add (1+ frame) = 0;
            end
        end
    end
end

nn = DATAPADDING_MSECS* (Fs/ 1000) + maxNsamples;
tweaked_deg = zeros( 1, nn);
% fprintf( 'nn is %d\n', nn);

for i= SEARCHBUFFER* Downsample+ 1: nn- SEARCHBUFFER* Downsample
    utt = Nutterances;

    while ((utt >= 1) && ((Utt_Start (utt)- 1)* Downsample> i))
        utt = utt- 1;
    end
    if (utt >= 1)
        delay = Utt_Delay (utt);
    else
        delay = Utt_Delay (1);
    end

    j = i + delay;
    if (j < SEARCHBUFFER * Downsample+ 1)
        j = SEARCHBUFFER * Downsample+ 1;
    end
    if (j > nn - SEARCHBUFFER * Downsample)
        j = nn - SEARCHBUFFER * Downsample;
    end
    tweaked_deg (i) = deg_data (j);
end

if (there_is_a_bad_frame)

    for frame = 0: stop_frame
        frame_is_bad (1+ frame) = (frame_disturbance (1+ frame)...
            > THRESHOLD_BAD_FRAMES);
        smeared_frame_is_bad (1+ frame) = FALSE;
    end
    frame_is_bad (1) = FALSE;
    SMEAR_RANGE = 2;

    for frame = SMEAR_RANGE: stop_frame- 1- SMEAR_RANGE
        max_itself_and_left = frame_is_bad (1+ frame);
        max_itself_and_right = frame_is_bad (1+ frame);

        for i = -SMEAR_RANGE: 0
            if (max_itself_and_left < frame_is_bad (1+ frame+ i))
                max_itself_and_left = frame_is_bad (1+ frame+ i);
            end
        end

        for i = 0: SMEAR_RANGE
            if (max_itself_and_right < frame_is_bad (1+ frame + i))
                max_itself_and_right = frame_is_bad (1+ frame + i);
            end
        end

        mini = max_itself_and_left;
        if (mini > max_itself_and_right)
            mini = max_itself_and_right;
        end

        smeared_frame_is_bad (1+ frame) = mini;
    end

    MINIMUM_NUMBER_OF_BAD_FRAMES_IN_BAD_INTERVAL = 5;
    number_of_bad_intervals = 0;
    frame = 0;
    while (frame <= stop_frame)
        while ((frame <= stop_frame) && (~smeared_frame_is_bad (1+ frame)))
            frame= frame+ 1;
        end

        if (frame <= stop_frame)
            start_frame_of_bad_interval(1+ number_of_bad_intervals)= ...
                1+ frame;

            while ((frame <= stop_frame) && (...
                    smeared_frame_is_bad (1+ frame)))
                frame= frame+ 1;
            end

            if (frame <= stop_frame)
                stop_frame_of_bad_interval(1+ number_of_bad_intervals)= ...
                    1+ frame;
                if (stop_frame_of_bad_interval(1+ number_of_bad_intervals)- ...
                        start_frame_of_bad_interval(1+ number_of_bad_intervals)...
                        >= MINIMUM_NUMBER_OF_BAD_FRAMES_IN_BAD_INTERVAL)
                    number_of_bad_intervals= number_of_bad_intervals+ 1;
                end
            end
        end
    end

    for bad_interval = 0: number_of_bad_intervals - 1
        start_sample_of_bad_interval(1+ bad_interval) = ...
            (start_frame_of_bad_interval(1+ bad_interval)- 1) * (Nf/ 2) ...
            + SEARCHBUFFER * Downsample+ 1;
        stop_sample_of_bad_interval(1+ bad_interval) = ...
            (stop_frame_of_bad_interval(1+ bad_interval)- 1) * (Nf/ 2) ...
            + Nf + SEARCHBUFFER* Downsample;
        if (stop_frame_of_bad_interval(1+ bad_interval) > stop_frame+ 1)
            stop_frame_of_bad_interval(1+ bad_interval) = stop_frame+ 1;
        end

        number_of_samples_in_bad_interval(1+ bad_interval) = ...
            stop_sample_of_bad_interval(1+ bad_interval) - ...
            start_sample_of_bad_interval(1+ bad_interval)+ 1;
    end
%     fprintf( 'number of bad intervals %d\n', number_of_bad_intervals);
%     fprintf( '%d %d\n', number_of_samples_in_bad_interval(1), ...
%         number_of_samples_in_bad_interval(2));
%     fprintf( '%d %d\n', start_sample_of_bad_interval(1), ...
%         start_sample_of_bad_interval(2));

    SEARCH_RANGE_IN_TRANSFORM_LENGTH = 4;
    search_range_in_samples= SEARCH_RANGE_IN_TRANSFORM_LENGTH * Nf;

    for bad_interval= 0: number_of_bad_intervals- 1
        ref = zeros (1, 2 * search_range_in_samples + ...
            number_of_samples_in_bad_interval (1+ bad_interval));
        deg = zeros (1, 2 * search_range_in_samples + ...
            number_of_samples_in_bad_interval (1+ bad_interval));

        ref(1: search_range_in_samples) = 0;

        ref (search_range_in_samples+ 1: search_range_in_samples+ ...
                number_of_samples_in_bad_interval (1+ bad_interval)) = ...
                ref_data (start_sample_of_bad_interval( 1+ bad_interval) + 1: ...
                start_sample_of_bad_interval( 1+ bad_interval) + ...
                number_of_samples_in_bad_interval (1+ bad_interval));

        ref (search_range_in_samples + ...
                number_of_samples_in_bad_interval (1+ bad_interval) + 1: ...
                search_range_in_samples + ...
                number_of_samples_in_bad_interval (1+ bad_interval) + ...
                search_range_in_samples) = 0;

        for i = 0: 2 * search_range_in_samples + ...
                number_of_samples_in_bad_interval (1+ bad_interval) - 1
            j = start_sample_of_bad_interval (1+ bad_interval) - ...
                search_range_in_samples + i;
            nn = maxNsamples - SEARCHBUFFER * Downsample + ...
                DATAPADDING_MSECS  * (Fs / 1000);
            if (j <= SEARCHBUFFER * Downsample)
                j = SEARCHBUFFER * Downsample+ 1;
            end
            if (j > nn)
                j = nn;
            end
            deg (1+ i) = tweaked_deg (j);
        end

        [delay_in_samples, best_correlation]= compute_delay ...
            (1, 2 * search_range_in_samples + ...
            number_of_samples_in_bad_interval (1+ bad_interval), ...
            search_range_in_samples, ref, deg);
        delay_in_samples_in_bad_interval (1+ bad_interval) =  ...
            delay_in_samples;
%         fprintf( 'delay_in_samples, best_correlation is \n\t%d, %f\n', ...
%             delay_in_samples, best_correlation);
%
        if (best_correlation < 0.5)
            delay_in_samples_in_bad_interval  (1+ bad_interval) = 0;
        end
    end

    if (number_of_bad_intervals > 0)
        doubly_tweaked_deg = tweaked_deg( 1: maxNsamples + ...
            DATAPADDING_MSECS  * (Fs / 1000));
        for bad_interval= 0: number_of_bad_intervals- 1
            delay = delay_in_samples_in_bad_interval (1+ bad_interval);

            for i = start_sample_of_bad_interval (1+ bad_interval): ...
                    stop_sample_of_bad_interval (1+ bad_interval)
                j = i + delay;
                if (j < 1)
                    j = 1;
                end
                if (j > maxNsamples)
                    j = maxNsamples;
                end
                h = tweaked_deg (j);
                doubly_tweaked_deg (i) = h;
            end
        end

        untweaked_deg = deg_data;
        deg_data = doubly_tweaked_deg;

        for bad_interval= 0: number_of_bad_intervals- 1
            for frame = start_frame_of_bad_interval (1+ bad_interval): ...
                    stop_frame_of_bad_interval (1+ bad_interval)- 1
                frame= frame- 1;
                start_sample_ref = SEARCHBUFFER * Downsample + ...
                    frame * Nf / 2+ 1;
                start_sample_deg = start_sample_ref;
                hz_spectrum_deg= short_term_fft (Nf, deg_data, ...
                    Whanning, start_sample_deg);
                pitch_pow_dens_deg( 1+ frame, :)= freq_warping (...
                    hz_spectrum_deg, Nb, frame);
            end

            oldScale = 1;
            for frame = start_frame_of_bad_interval (1+ bad_interval): ...
                    stop_frame_of_bad_interval (1+ bad_interval)- 1
                frame= frame- 1;
                % see implementation for detail why 1 needed to be
                % subtracted
                total_audible_pow_ref = total_audible (frame, ...
                    pitch_pow_dens_ref, 1);
                total_audible_pow_deg = total_audible (frame, ...
                    pitch_pow_dens_deg, 1);
                scale = (total_audible_pow_ref + 5e3) / ...
                    (total_audible_pow_deg + 5e3);
                if (frame > 0)
                    scale = 0.2 * oldScale + 0.8*scale;
                end
                oldScale = scale;
                if (scale > MAX_SCALE)
                    scale = MAX_SCALE;
                end
                if (scale < MIN_SCALE)
                    scale = MIN_SCALE;
                end

                pitch_pow_dens_deg (1+ frame, :) = ...
                    pitch_pow_dens_deg (1+ frame, :)* scale;
                loudness_dens_ref= intensity_warping_of (frame, ...
                    pitch_pow_dens_ref);
                loudness_dens_deg= intensity_warping_of (frame, ...
                    pitch_pow_dens_deg);
                disturbance_dens = loudness_dens_deg - loudness_dens_ref;

                for band = 1: Nb
                    deadzone(band) = min (loudness_dens_deg(band), ...
                        loudness_dens_ref(band));
                    deadzone(band) = deadzone(band)* 0.25;
                end

                for band = 1: Nb
                    d = disturbance_dens (band);
                    m = deadzone (band);

                    if (d > m)
                        disturbance_dens (band) = ...
                            disturbance_dens (band)- m;
                    else
                        if (d < -m)
                            disturbance_dens (band) = ...
                                disturbance_dens (band)+ m;
                        else
                            disturbance_dens (band) = 0;
                        end
                    end
                end

                frame_disturbance( 1+ frame) = min (...
                    frame_disturbance( 1+ frame), pseudo_Lp(...
                    disturbance_dens, D_POW_F));
                disturbance_dens= multiply_with_asymmetry_factor ...
                    (disturbance_dens, frame, pitch_pow_dens_ref, ...
                    pitch_pow_dens_deg);
                frame_disturbance_asym_add(1+ frame) = min (...
                    frame_disturbance_asym_add(1+ frame), ...
                    pseudo_Lp (disturbance_dens, A_POW_F));
            end
        end
        deg_data = untweaked_deg;
    end
end

for frame = 0: stop_frame
    h = 1;
    if (stop_frame + 1 > 1000)
        n = floor( (maxNsamples - 2 * SEARCHBUFFER * Downsample)...
            / (Nf / 2)) - 1;
        timeWeightFactor = (n - 1000) / 5500;
        if (timeWeightFactor > 0.5)
            timeWeightFactor = 0.5;
        end
        h = (1.0 - timeWeightFactor) + timeWeightFactor * frame / n;
    end

    time_weight (1 +frame) = h;
end

% fid= fopen( 'tmp_mat1.txt', 'at');
% fprintf( '\n');
for frame = 0: stop_frame
    h = ((total_power_ref (1+ frame) + 1e5) / 1e7)^ 0.04;
%     if (frame== 118)
%         fprintf( '%f\n', h);
%         fprintf( '%f\n', frame_disturbance( 1+ frame));
%     end
    frame_disturbance( 1+ frame) = frame_disturbance( 1+ frame)/ h;

%     if (frame== 118)
%         fprintf( '%f\n', frame_disturbance( 1+ frame));
%     end
%
    frame_disturbance_asym_add( 1+ frame) = ...
        frame_disturbance_asym_add( 1+ frame)/ h;
    if (frame_disturbance( 1+ frame) > 45)
        frame_disturbance( 1+ frame) = 45;
    end
    if (frame_disturbance_asym_add( 1+ frame)> 45)
        frame_disturbance_asym_add( 1+ frame) = 45;
    end
end
% fclose ( fid);

d_indicator = Lpq_weight (start_frame, stop_frame, ...
    D_POW_S, D_POW_T, frame_disturbance, time_weight);
a_indicator = Lpq_weight (start_frame, stop_frame, ...
    A_POW_S, A_POW_T, frame_disturbance_asym_add, time_weight);

pesq_mos = 4.5 - D_WEIGHT * d_indicator - A_WEIGHT * a_indicator;

if (Plot_Frame> 0)
    figure;
    subplot( 1, 2, 1);
    mesh( 0: stop_frame, centre_of_band_hz, D_disturbance');
    title( 'disturbance');
    subplot( 1, 2, 2);
    mesh( 0: stop_frame, centre_of_band_hz, DA_disturbance');
    title( 'disturbance after asymmetry processing');
end

% fid= fopen( 'tmp_mat.txt', 'wt');
% fprintf( fid, 'time_weight\n');
% fprintf( fid, '%f\n', time_weight);
% fprintf( fid, 'frame_disturbance:\n');
% fprintf( fid, '%f\n', frame_disturbance);
% fprintf( fid, 'frame_disturbance_asym_add\n');
% fprintf( fid, '%f\n', frame_disturbance_asym_add);
% fclose( fid);

function result_time= Lpq_weight(start_frame, stop_frame, ...
        power_syllable, power_time, frame_disturbance, time_weight)

global NUMBER_OF_PSQM_FRAMES_PER_SYLLABE

% fid= fopen( 'tmp_mat1.txt', 'at');
% fprintf( 'result_time:\n');

result_time= 0;
total_time_weight_time = 0;
% fprintf( 'start/end frame: %d/%d\n', start_frame, stop_frame);
for start_frame_of_syllable = start_frame: ...
        NUMBER_OF_PSQM_FRAMES_PER_SYLLABE/2: stop_frame
    result_syllable = 0;
    count_syllable = 0;

    for frame = start_frame_of_syllable: ...
            start_frame_of_syllable + NUMBER_OF_PSQM_FRAMES_PER_SYLLABE- 1
        if (frame <= stop_frame)
            h = frame_disturbance(1+ frame);
%             if (start_frame_of_syllable== 101)
%                 fprintf( fid, '%f\n', h);
%             end
            result_syllable = result_syllable+ (h^ power_syllable);
        end
        count_syllable = count_syllable+ 1;
    end

    result_syllable = result_syllable/ count_syllable;
    result_syllable = result_syllable^ (1/power_syllable);

    result_time= result_time+ (time_weight (...
        1+ start_frame_of_syllable - start_frame) * ...
        result_syllable)^ power_time;
    total_time_weight_time = total_time_weight_time+ ...
        time_weight (1+ start_frame_of_syllable - start_frame)^ power_time;

%     fprintf( fid, '%f\n', result_time);
end
% fclose (fid);

% fprintf( 'total_time_weight_time is %f\n', total_time_weight_time);
result_time = result_time/ total_time_weight_time;
result_time= result_time^ (1/ power_time);
% fprintf( 'result_time is %f\n\n', result_time);

function [best_delay, max_correlation] = compute_delay (...
    start_sample, stop_sample, search_range, ...
    time_series1, time_series2)

n = stop_sample - start_sample+ 1;
power_of_2 = 2^ (ceil( log2( 2 * n)));

power1 = pow_of (time_series1, start_sample, stop_sample, n)* ...
    n/ power_of_2;
power2 = pow_of (time_series2, start_sample, stop_sample, n)* ...
    n/ power_of_2;
normalization = sqrt (power1 * power2);
% fprintf( 'normalization is %f\n', normalization);

if ((power1 <= 1e-6) || (power2 <= 1e-6))
    max_correlation = 0;
    best_delay= 0;
end

x1( 1: power_of_2)= 0;
x2( 1: power_of_2)= 0;
y( 1: power_of_2)= 0;

x1( 1: n)= abs( time_series1( start_sample: ...
    stop_sample));
x2( 1: n)= abs( time_series2( start_sample: ...
    stop_sample));

x1_fft= fft( x1, power_of_2)/ power_of_2;
x2_fft= fft( x2, power_of_2);
x1_fft_conj= conj( x1_fft);
y= ifft( x1_fft_conj.* x2_fft, power_of_2);

best_delay = 0;
max_correlation = 0;

% these loop can be rewritten
for i = -search_range: -1
    h = abs (y (1+ i + power_of_2)) / normalization;
    if (h > max_correlation)
        max_correlation = h;
        best_delay= i;
    end
end
for i = 0: search_range- 1
    h = abs (y (1+i)) / normalization;
    if (h > max_correlation)
        max_correlation = h;
        best_delay= i;
    end
end
best_delay= best_delay- 1;

function mod_disturbance_dens= multiply_with_asymmetry_factor (...
    disturbance_dens, frame, pitch_pow_dens_ref, pitch_pow_dens_deg)

global Nb
for i = 1: Nb
    ratio = (pitch_pow_dens_deg(1+ frame, i) + 50)...
        / (pitch_pow_dens_ref (1+ frame, i) + 50);
    h = ratio^ 1.2;
    if (h > 12)
        h = 12;
    elseif (h < 3)
        h = 0.0;
    end
    mod_disturbance_dens (i) = disturbance_dens (i) * h;
end

function loudness_dens = intensity_warping_of (...
    frame, pitch_pow_dens)

global abs_thresh_power Sl Nb centre_of_band_bark
ZWICKER_POWER= 0.23;
for band = 1: Nb
    threshold = abs_thresh_power (band);
    input = pitch_pow_dens (1+ frame, band);

    if (centre_of_band_bark (band) < 4)
        h =  6 / (centre_of_band_bark (band) + 2);
    else
        h = 1;
    end

    if (h > 2)
        h = 2;
    end
    h = h^ 0.15;
    modified_zwicker_power = ZWICKER_POWER * h;
    if (input > threshold)
        loudness_dens (band) = ((threshold / 0.5)^ modified_zwicker_power)...
            * ((0.5 + 0.5 * input / threshold)^ modified_zwicker_power- 1);
    else
        loudness_dens (band) = 0;
    end

    loudness_dens (band) = loudness_dens (band)* Sl;
end

function result= pseudo_Lp (x, p)

global Nb width_of_band_bark
totalWeight = 0;
result = 0;
for band = 2: Nb
    h = abs (x (band));
    w = width_of_band_bark (band);
    prod = h * w;

    result = result+ prod^ p;
    totalWeight = totalWeight+ w;
end
result = (result/ totalWeight)^ (1/p);
result = result* totalWeight;

function mod_pitch_pow_dens_ref= freq_resp_compensation (number_of_frames, ...
    pitch_pow_dens_ref, avg_pitch_pow_dens_ref, ...
    avg_pitch_pow_dens_deg, constant)

global Nb

for band = 1: Nb
    x = (avg_pitch_pow_dens_deg (band) + constant) / ...
        (avg_pitch_pow_dens_ref (band) + constant);
    if (x > 100.0)
        x = 100.0;
    elseif (x < 0.01)
        x = 0.01;
    end

    for frame = 1: number_of_frames
        mod_pitch_pow_dens_ref(frame, band) = ...
            pitch_pow_dens_ref(frame, band) * x;
    end
end

function avg_pitch_pow_dens= time_avg_audible_of(number_of_frames, ...
    silent, pitch_pow_dens, total_number_of_frames)

global Nb abs_thresh_power

for band = 1: Nb
    result = 0;
    for frame = 1: number_of_frames
        if (~silent (frame))
            h = pitch_pow_dens (frame, band);
            if (h > 100 * abs_thresh_power (band))
                result = result + h;
            end
        end

        avg_pitch_pow_dens (band) = result/ total_number_of_frames;
    end
end

function hz_spectrum= short_term_fft (Nf, data, Whanning, start_sample)

x1= data( start_sample: start_sample+ Nf-1).* Whanning;
x1_fft= fft( x1);
hz_spectrum= abs( x1_fft( 1: Nf/ 2)).^ 2;
hz_spectrum( 1)= 0;

function pitch_pow_dens= freq_warping( hz_spectrum, Nb, frame)

global nr_of_hz_bands_per_bark_band pow_dens_correction_factor
global Sp

hz_band = 1;
for bark_band = 1: Nb
    n = nr_of_hz_bands_per_bark_band (bark_band);
    sum = 0;
    for i = 1: n
        sum = sum+ hz_spectrum( hz_band);
        hz_band= hz_band+ 1;
    end
    sum = sum* pow_dens_correction_factor (bark_band);
    sum = sum* Sp;
    pitch_pow_dens (bark_band) = sum;

end

function total_audible_pow = total_audible (frame, ...
    pitch_pow_dens, factor)

global Nb abs_thresh_power

total_audible_pow = 0;
for band= 2: Nb
    h = pitch_pow_dens (frame+ 1,band);
    threshold = factor * abs_thresh_power (band);
    if (h > threshold)
        total_audible_pow = total_audible_pow+ h;
    end
end

author: vipchengrui

repo: traditional-speech-enhancement


text:
%  Copyright (C) 2004 Yokogawa Electric Corporation, All Rights Reserved.
%
%  usage:
%    calib_info = GetMeg160CalibInfoM( fid )
%
%  arguments:
%    fid             : file ID
%
%  return values:
%    calib_info      : (m x 3) matrix of calibration information. (m: channel_count)
%                        detail of n: (If the sensor doesn't have following item, the item's entry will be 'Inf'.)
%                            1   : channel number in Meg160(zero start)
%                            2   : gain [Tesla/V]
%                            3   : offset voltage [V]
%
%  confirmation of revision:
%   GetMeg160CalibInfoM( Inf ) will show and return revision of this function.
%

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [success path] = traverse_world(start,pot)
%
% traverse world
%
res(size(pot,1),size(pot,2))=0;
t = 1;
nr = start.r;
nc = start.c;
res(nr,nc,t)=1;
while((t<size(pot,3))&&(pot(nr,nc,t)~=0)&&(pot(nr,nc,t)~=1))
    if(nr ~= 1)
        n = pot(nr-1,nc,t+1);
    else
        n = Inf;
    end
    if(nc ~= size(pot,2))
        e = pot(nr,nc+1,t+1);
    else
        e = Inf;
    end
    if(nr ~= size(pot,1))
        s = pot(nr+1,nc,t+1);
    else
        s = Inf;
    end
    if(nc ~= 1)
        w = pot(nr,nc-1,t+1);
    else
        w = Inf;
    end
    st = pot(nr,nc,t+1);
    minv = min([n,e,s,w,st]);
    if (minv == st)
        %do nothing
    elseif((minv == n)&&(nr ~= 1))
        nr = nr - 1;
    elseif((minv == e)&&(nc ~= 10))
        nc = nc + 1;
    elseif((minv == s)&&(nr ~= 10))
        nr = nr + 1;
    elseif((minv == w)&&(nc ~= 1))
        nc = nc - 1;
    end
    t = t + 1;
    res(nr,nc,t)=1;
end
if(pot(nr,nc,t)==0)
    success = 1; %if goal reached
else
    success = 0; %if goal no longer reachable
end
path = res;
end

author: Sable

repo: mcbench-benchmarks


text:
function [ readerobj ] = subset_reader(ro, dim1_limits, dim2_limits)
%SUBSET_READER Creates a reader with visibility into only part of a file
%
% Written by: Wade Schwartzkopf, NGA Research
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

%% Define object methods
readerobj = ro; % Copy most methods
readerobj.read_chip = @chip_fun; % Subset chipper
% Subset metadata
meta = readerobj.get_meta();
meta.ImageData.NumCols = diff(dim1_limits) + 1;
meta.ImageData.NumRows = diff(dim2_limits) + 1;
if isfield(meta.ImageData,'FirstCol')
    meta.ImageData.FirstCol = meta.ImageData.FirstCol + dim1_limits(1) - 1;
end
if isfield(meta.ImageData,'FirstRow')
    meta.ImageData.FirstRow = meta.ImageData.FirstRow + dim2_limits(1) - 1;
end
readerobj.get_meta = @() meta;

    % Check to see if indices are out of bounds of subset area and index
    % based on subset, not full image.
    function out = chip_fun(varargin)
        [ dim1range, dim2range, subsample ] =...
            check_chipper_args( [diff(dim1_limits) diff(dim2_limits)] + 1, varargin{:});
        out = ro.read_chip(dim1range + double(dim1_limits(1)) - 1, ...
            dim2range + double(dim2_limits(1)) - 1, subsample, varargin{4:end});
    end

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
% Small technical example.
% Shows how to retrieve the details on a rectangular grid, after a multi-
% resolution decomposition has been performed.
%
disp('Small technical example.');
disp('Shows how to retrieve the details on a rectangular grid, after a multi-');
disp('resolution decomposition has been performed.');
disp('FOR MORE INFORMATION:  help retrieveR');
disp(' ');
disp('See also the report http://repository.cwi.nl:8888/cwi_repository/docs/IV/04/04178D.pdf');
disp('Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>');
disp(' (C) 1998-2006 Stichting CWI, Amsterdam, The Netherlands');
disp(' ');
%---PARAMETERS-----------------------------------------------------------------
% How to execute, set parameters
N = 6;                   %  maximum level (even number) in lifting scheme
filtername = 'Neville4';
%
%---INSERT YOUR IMAGE HERE-----------------------------------------------------
if exist('imread','file') == 2
  Orig = double(imread('zenithgray.TIF','tiff'));
else
  load zenithgray; Orig = zenithgray; clear zenithgray;
end
%
disp([' Dimensions of original      ' int2str( size(Orig) )]);
%---DECOMPOSITION--------------------------------------------------------------
disp([' Filter type is ' filtername]);
[C,S] = QLiftDec2(Orig,N,filtername);
%
%---RETRIEVE APPROXIMATION-----------------------------------------------------
A = retrieveR(N, 'a', C, S);
sizeA = size(A);
%
%---RETRIEVE DETAIL AT EVEN LEVEL----------------------------------------------
level = N-2;
disp([' Level ' int2str(level)]);
D = retrieveR(level, 'd', C, S);
%
disp(['The dimensions of the detail function read as follows ' int2str(size(D))]);


text:
function tstoolInstall
% Install OpenTSTOOL by saving it in Matlab search path

if (  isempty(which('units.mat')) ...
   || isempty(which('nn_prepare')) ...
   || isempty(which('genbyode')) ...
   || isempty(which('opentstool')) )
	error('Seems like OpenTSTOOL directories are not in current search path. Run tstoolInit first')
else
    if ( usejava('awt') )
	h=msgbox(['I will start the Matlab path tool. Just click `save` there to save the current Matlab path.' ...
	'If you do not have permission to save the global Matlab search path, save it ' ...
	'in your userpath ' userpath ' instead.'], 'OpenTSTOOL installation','warn','modal' );
        uiwait(h);
    else
        fprintf(['I will now start the path tool. Just click `save` to save the current Matlab path.\n' ...
	'If you do not have permission to save the global Matlab search path, save it\n' ...
	'in your userpath %s instead\n\nPress a key to continue...\n'],userpath)
	pause
    end
    pathtool
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function res = isfield(this, varargin)
% Returns true if the string fieldname is the name of a field in the
% substructure 'other' in the meeg object 'this'.
% FORMAT res = isfield(this,fieldname)
%
% An overloaded function...
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Christophe Phillips
% $Id: isfield.m 2720 2009-02-09 19:50:46Z vladimir $

res = isfield(this.other, varargin{:});


text:
function [x,y] = EquiNodes2D(N);

% function [x,y] = EquiNodes2D(N);
% Purpose  : Compute (x,y) nodes in equilateral triangle for polynomial of order N

% total number of nodes
Np = (N+1)*(N+2)/2;

% Create equidistributed nodes on equilateral triangle
L1 = zeros(Np,1); L2 = zeros(Np,1); L3 = zeros(Np,1);
sk = 1;
for n=1:N+1
  for m=1:N+2-n
    L1(sk) = (n-1)/N; L3(sk) = (m-1)/N;
    sk = sk+1;
  end
end
L2 = 1.0-L1-L3;

x = -L2+L3; y = (-L2-L3+2*L1)/sqrt(3.0);
return

repo: Aero-matlab


text:
% File: FitG.m
%
% Copyright (C) Daphne Koller, Stanford Univerity, 2012

function [mu sigma] = FitG(X, W)

% X: (N x 1): N examples (1 dimensional)
% W: (N x 1): Weights over examples (W(i) is the weight for X(i))

% Fit N(mu, sigma^2) to the empirical distribution

mu = 0;
sigma = 1;

mu = W'*X/sum(W);
v = W'*(X.*X)/sum(W) - mu^2;
sigma = sqrt(v);

author: anhncs

repo: Probabilistic-Graphical-Models


text:
%RODRIGUES  Converts a rotation matrix to a rotation vector or vice versa
%
%     dst = cv.Rodrigues(src)
%     [dst,jacobian] = cv.Rodrigues(src)
%
% ## Input
% * __src__ Input rotation vector (3x1 or 1x3) or rotation matrix (3x3). Both
%   single and double-precision floating-point types are supported.
%
% ## Output
% * __dst__ Output rotation matrix (3x3) or rotation vector (3x1 or 1x3),
%   respectively. Same data type as `src`.
% * __jacobian__ Optional output Jacobian matrix, 3x9 or 9x3, which is a
%   matrix of partial derivatives of the output array components with respect
%   to the input array components. Same data type as `src`.
%
% The function transforms a rotation matrix in the following way:
%
%     theta <- norm(r)
%     r <- r/theta
%     R = cos(theta) * I + (1 - cos(theta)) * r * r^T + sin(theta) * A
%     A = [0, -rz, ry; rz, 0, -rx; -ry, rx, 0]
%
% Inverse transformation can be also done easily, since
%
%     sin(theta) * A = (R - R^T) / 2
%
% A rotation vector is a convenient and most compact representation of a
% rotation matrix (since any rotation matrix has just 3 degrees of
% freedom). The representation is used in the global 3D geometry
% optimization procedures like cv.calibrateCamera, cv.stereoCalibrate, or
% cv.solvePnP.
%
% See also: cv.calibrateCamera, cv.stereoCalibrate, cv.solvePnP,
%  rotationMatrixToVector, rotationVectorToMatrix
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [Q, V, policy, mean_discrepancy] = mdp_Q_learning(P, R, discount, N)

% mdp_Q_learning   Evaluation of the matrix Q, using the Q learning algorithm
%
% Arguments
% -------------------------------------------------------------------------
% Let S = number of states, A = number of actions
%   P(SxSxA)  = transition matrix
%              P could be an array with 3 dimensions or
%              a cell array (1xA), each cell containing a sparse matrix (SxS)
%   R(SxSxA) or (SxA) = reward matrix
%              R could be an array with 3 dimensions (SxSxA) or
%              a cell array (1xA), each cell containing a sparse matrix (SxS) or
%              a 2D array(SxA) possibly sparse
%   discount  = discount rate in ]0; 1[
%   N(optional) = number of iterations to execute, default value: 10000.
%                 It is an integer greater than the default value.
% Evaluation --------------------------------------------------------------
%   Q(SxA) = learned Q matrix
%   V(S)   = learned value function.
%   policy(S) = learned optimal policy.
%   mean_discrepancy(N/100) = vector of V discrepancy mean over 100 iterations
%             Then the length of this vector for the default value of N is 100.

% MDPtoolbox: Markov Decision Processes Toolbox
% Copyright (C) 2009  INRA
% Redistribution and use in source and binary forms, with or without modification,
% are permitted provided that the following conditions are met:
%    * Redistributions of source code must retain the above copyright notice,
%      this list of conditions and the following disclaimer.
%    * Redistributions in binary form must reproduce the above copyright notice,
%      this list of conditions and the following disclaimer in the documentation
%      and/or other materials provided with the distribution.
%    * Neither the name of the <ORGANIZATION> nor the names of its contributors
%      may be used to endorse or promote products derived from this software
%      without specific prior written permission.
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
% ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
% WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
% IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
% INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
% BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
% DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
% OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
% OF THE POSSIBILITY OF SUCH DAMAGE.

% check of arguments
if (discount <= 0 || discount >= 1)
    disp('--------------------------------------------------------')
    disp('MDP Toolbox ERROR: Discount rate must be in ]0,1[')
    disp('--------------------------------------------------------')
elseif (nargin >= 4) && (N < 10000)
    disp('--------------------------------------------------------')
    disp('MDP Toolbox ERROR: N must be upper than 10000')
    disp('--------------------------------------------------------')
else

    % initialization of optional arguments
    if (nargin < 4); N=10000; end;

    % Find number of states and actions
    if iscell(P)
        S = size(P{1},1);
        A = length(P);
    else
        S = size(P,1);
        A = size(P,3);
    end;

    % Initialisations
    Q = zeros(S,A);
    dQ = zeros(S,A);
    mean_discrepancy = [];
    discrepancy = [];

    % Initial state choice
    s = randi([1,S]);

    for n=1:N

        % Reinitialisation of trajectories every 100 transitions
        if (mod(n,100)==0); s = randi([1,S]); end;

        % Action choice : greedy with increasing probability
        % probability 1-(1/log(n+2)) can be changed
        pn = rand(1);
        if (pn < (1-(1/log(n+2))))
          [nil,a] = max(Q(s,:));
        else
          a = randi([1,A]);
        end;

        % Simulating next state s_new and reward associated to <s,s_new,a>
        p_s_new = rand(1);
        p = 0;
        s_new = 0;
        while ((p < p_s_new) && (s_new < S))
            s_new = s_new+1;
            if iscell(P)
                p = p + P{a}(s,s_new);
            else
                p = p + P(s,s_new,a);
            end;
        end;
        if iscell(R)
            r = R{a}(s,s_new);
        elseif ndims(R) == 3
            r = R(s,s_new,a);
        else
            r = R(s,a);
        end;

        % Updating the value of Q
        % Decaying update coefficient (1/sqrt(n+2)) can be changed
        delta = r + discount*max(Q(s_new,:)) - Q(s,a);
        dQ = (1/sqrt(n+2))*delta;
        Q(s,a) = Q(s,a) + dQ;

        % Current state is updated
        s = s_new;

        % Computing and saving maximal values of the Q variation
        discrepancy(mod(n,100)+1) = abs(dQ);

        % Computing means all over maximal Q variations values
        if (length(discrepancy) == 100)
           mean_discrepancy = [ mean_discrepancy mean(discrepancy)];
           discrepancy = [];
        end;

    end;

    %compute the value function and the policy
    [V, policy] = max(Q,[],2);

end;


text:
function filterSteerable( theta )
% Steerable 2D Gaussian derivative filter (for visualization).
%
% This function is a demonstration of steerable filters.  The directional
% derivative of G in an arbitrary direction theta can be found by taking a
% linear combination of the directional derivatives dxG and dyG.
%
% USAGE
%  filterSteerable( theta )
%
% INPUTS
%  theta   - orientation in radians
%
% OUTPUTS
%
% EXAMPLE
%  filterSteerable( pi/4 );
%
% See also filterGauss
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.0
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

% Get G
[x,y]=meshgrid(-1:.1:1, -1:.1:1 );
r = sqrt( x.^2 + y.^2 );
G = exp( -r .* r *2 );

% get first derivatives of G.  note: d/dx(G)=-2x*G
phi = atan2( y, x );
dxG = r .* cos(phi) .* G;
dyG = r .* sin(phi) .* G;

% get directional derivative by taking linear comb in theta
Gtheta = cos(theta)*dxG + sin(theta)*dyG;

% dislpay (scale for visualization purposes)
GS = cat(3,G,dxG*2,dyG*2,Gtheta*2);
figure(1); montage2( GS );

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function [mustUL, pos_mustUL, mustUL_linear, pos_mustUL_linear] = findMustULWithGAMS(model, minFluxesW, maxFluxesW, constrOpt, excludedRxns, mustSetFirstOrder, solverName, runID, outputFolder, outputFileName, printExcel, printText, printReport, keepInputs, keepGamsOutputs, verbose)
% This function runs the second step of `optForce`, that is to solve a
% bilevel mixed integer linear programming  problem to find a second order
% MustUL set.
%
% USAGE:
%
%    [mustUL, pos_mustUL, mustUL_linear, pos_mustUL_linear] = findMustULWithGAMS(model, minFluxesW, maxFluxesW, varargin)
%
% INPUTS:
%    model:                      (structure) a metabolic model with at
%                                least the following fields:
%
%                                  * .rxns - Reaction IDs in the model
%                                  * .mets - Metabolite IDs in the model
%                                  * .S - Stoichiometric matrix (sparse)
%                                  * .b - RHS of `Sv = b` (usually zeros)
%                                  * .c - Objective coefficients
%                                  * .lb - Lower bounds for fluxes
%                                  * .ub - Upper bounds for fluxes
%    minFluxesW:                 (double array of size `n_rxns x 1`) minimum
%                                fluxes for each reaction in the model for
%                                wild-type strain. This can be obtained by
%                                running the function `FVAOptForce`. E.g.:
%                                `minFluxesW = [-90; -56];`
%    maxFluxesW:                 (double array of size n_rxns x 1) maximum
%                                fluxes for each reaction in the model for
%                                wild-type strain. This can be obtained by
%                                running the function `FVAOptForce`. E.g.:
%                                `maxFluxesW = [90; 56];`
%
% OPTIONAL INPUTS:
%    constrOpt:                  (Structure) structure containing
%                                additional contraints. Include here only
%                                reactions whose flux is fixed, i.e.,
%                                reactions whose lower and upper bounds
%                                have the same value. Do not include here
%                                reactions whose lower and upper bounds
%                                have different values. Such contraints
%                                should be defined in the lower and upper
%                                bounds of the model. The structure has the
%                                following fields:
%
%                                  * .rxnList - Reaction list (cell array)
%                                  * .values -  Values for constrained
%                                    reactions (double array)
%                                    E.g.: `struct('rxnList', {{'EX_gluc', 'R75', 'EX_suc'}}, 'values', [-100, 0, 155.5]');`
%    excludedRxns:               (cell array) Reactions to be excluded to
%                                the `MustUL` set. This could be used to
%                                avoid finding transporters or exchange
%                                reactions in the set. Default = empty.
%    mustSetFirstOrder:          (cell array) Reactions that belong to
%                                `MustU` and `MustL` (first order sets).
%                                Default = empty.
%    solverName:                 (string) Name of the solver used in
%                                GAMS. Default = 'cplex'.
%    runID:                      (string) ID for identifying this run.
%                                Default = ['run' date hour].
%    outputFolder:               (string) name for folder in which
%                                results will be stored.
%                                Default = 'OutputsFindMustUL'.
%    outputFileName:             (string) name for files in which
%                                results. will be stored
%                                Default = 'MustULSet'.
%    printExcel:                 (double) boolean to describe wheter
%                                data must be printed in an excel file or
%                                not. Default = 1
%    printText:                  (double) boolean to describe wheter
%                                data must be printed in an plaint text
%                                file or not. Default = 1
%    printReport:                (double) 1 to generate a report in a
%                                plain text file. 0 otherwise. Default = 1
%    keepInputs:                 (double) 1 to mantain folder with
%                                inputs to run `findMustUL.gms`. 0 otherwise.
%                                Default = 1
%    keepGamsOutputs:            (double) 1 to mantain files returned by
%                                `findMustUL.gms`. 0 otherwise. Default = 1
%    verbose:                    (double) 1 to print results in console.
%                                0 otherwise. Default = 0
%
% OUTPUTS:
%    mustUL:                     (cell array of size number of sets found X
%                                2) Cell array containing the reactions IDs
%                                which belong to the `MustUL` set. Each row
%                                contain a couple of reactions that must
%                                decrease their flux.
%    pos_mustUL:                 (double array of size number of sets found
%                                X 2) double array containing the positions
%                                of each reaction in `mustUL` with regard to
%                                model.rxns
%    mustUL_linear:              (cell array of size number of unique
%                                reactions found X 1) Cell array containing
%                                the unique reactions ID which belong to
%                                the `MustUL` Set
%    pos_mustUL_linear:          (double array of size number of unique
%                                reactions found X 1) double array
%                                containing positions for reactions in
%                                mustUL_linear. with regard to `model.rxns`
%    outputFileName.xls:         (file) File containing one column
%                                array with identifiers for reactions in
%                                MustUL. This file will only be generated
%                                if the user entered `printExcel = 1`. Note
%                                that the user can choose the name of this
%                                file entering the input `outputFileName = 'PutYourOwnFileNameHere';`
%    outputFileName.txt:         (file) File containing one column
%                                array with identifiers for reactions in
%                                MustUL. This file will only be generated
%                                if the user entered `printText = 1`. Note
%                                that the user can choose the name of this
%                                file entering the input `outputFileName = 'PutYourOwnFileNameHere';`
%    outputFileName_Info.xls:    (file) File containing one column
%                                array. In each row the user will find a
%                                couple of reactions. Each couple of
%                                reaction was found in one iteration of
%                                `FindMustUL.gms`. This file will only be
%                                generated if the user entered `printExcel = 1`.
%                                Note that the user can choose the name
%                                of this file entering the input
%                                `outputFileName = 'PutYourOwnFileNameHere';`
%    outputFileName_Info.txt:    (file) File containing one column
%                                array. In each row the user will find a
%                                couple of reactions. Each couple of
%                                reaction was found in one iteration of
%                                `FindMustUL.gms`. This file will only be
%                                generated if the user entered `printText = 1`.
%                                Note that the user can choose the name
%                                of this file entering the input
%                                `outputFileName = 'PutYourOwnFileNameHere';`
%    findMustUL.lst:             (file) file autogenerated by GAMS. It
%                                contains information about equations,
%                                variables, parameters as well as
%                                information about the running (values at
%                                each iteration). This file only will be
%                                saved in the output folder is the user
%                                entered `keepGamsOutputs = 1`
%    GtoMUL.gdx:                 (file) file containing values for
%                                variables, parameters, etc. which were
%                                found by GAMS when solving `findMustUL.gms`.
%                                This file only will be saved in the output
%                                folder is the user entered `keepInputs = 1`
%
% NOTE:
%
%    This function is based in the GAMS files written by Sridhar
%    Ranganathan which were provided by the research group of Costas D.
%    Maranas. For a detailed description of the `optForce` procedure, please
%    see: `Ranganathan S, Suthers PF, Maranas CD (2010) OptForce: An
%    Optimization Procedure for Identifying All Genetic Manipulations
%    Leading to Targeted Overproductions. PLOS Computational Biology 6(4):
%    e1000744`. https://doi.org/10.1371/journal.pcbi.1000744
%
% .. Author: - Sebastian Mendoza, May 30th 2017, Center for Mathematical Modeling, University of Chile, snmendoz@uc.cl

optionalParameters = {'constrOpt', 'excludedRxns', 'mustSetFirstOrder', 'solverName', 'runID', 'outputFolder', 'outputFileName',  ...
    'printExcel', 'printText', 'printReport', 'keepInputs', 'keepGamsOutputs', 'verbose'};

if (numel(varargin) > 0 && (~ischar(varargin{1}) || ~any(ismember(varargin{1},optionalParameters))))

    tempargin = cell(1,2*(numel(varargin)));
    for i = 1:numel(varargin)

        tempargin{2*(i-1)+1} = optionalParameters{i};
        tempargin{2*(i-1)+2} = varargin{i};
    end
    varargin = tempargin;

end

parser = inputParser();
parser.addRequired('model', @(x) isstruct(x) && isfield(x, 'S') && isfield(model, 'rxns')...
    && isfield(model, 'mets') && isfield(model, 'lb') && isfield(model, 'ub') && isfield(model, 'b')...
    && isfield(model, 'c'))
parser.addRequired('minFluxesW', @isnumeric)
parser.addRequired('maxFluxesW', @isnumeric)
parser.addParamValue('constrOpt', struct('rxnList', {{}}, 'values', []),@ (x) isstruct(x) && isfield(x, 'rxnList') && isfield(x, 'values') ...
    && length(x.rxnList) == length(x.values) && length(intersect(x.rxnList, model.rxns)) == length(x.rxnList))
parser.addParamValue('excludedRxns', {}, @(x) iscell(x) && length(intersect(x, model.rxns)) == length(x))
parser.addParamValue('mustSetFirstOrder', {}, @(x) iscell(x) && length(intersect(x, model.rxns)) == length(x))
solvers = checkGAMSSolvers('MIP');
if isempty(solvers)
    error('there is no GAMS solvers available to solve Mixed Integer Programming problems') ;
else
    if ismember('cplex', lower(solvers))
        defaultSolverName = 'cplex';
    else
        defaultSolverName = lower(solvers(1));
    end
end

parser.addParamValue('solverName', defaultSolverName, @(x) ischar(x))
hour = clock; defaultRunID = ['run-' date '-' num2str(hour(4)) 'h' '-' num2str(hour(5)) 'm'];
parser.addParamValue('runID', defaultRunID, @(x) ischar(x))
parser.addParamValue('outputFolder', 'OutputsFindMustUL', @(x) ischar(x))
parser.addParamValue('outputFileName', 'MustULSet', @(x) ischar(x))
parser.addParamValue('printExcel', 1, @(x) isnumeric(x) || islogical(x));
parser.addParamValue('printText', 1, @(x) isnumeric(x) || islogical(x));
parser.addParamValue('printReport', 1, @(x) isnumeric(x) || islogical(x));
parser.addParamValue('keepInputs', 1, @(x) isnumeric(x) || islogical(x));
parser.addParamValue('keepGamsOutputs', 1, @(x) isnumeric(x) || islogical(x));
parser.addParamValue('verbose', 1, @(x) isnumeric(x) || islogical(x));

parser.parse(model, minFluxesW, maxFluxesW, varargin{:})
model = parser.Results.model;
minFluxesW = parser.Results.minFluxesW;
maxFluxesW = parser.Results.maxFluxesW;
constrOpt= parser.Results.constrOpt;
excludedRxns= parser.Results.excludedRxns;
mustSetFirstOrder = parser.Results.mustSetFirstOrder;
solverName = parser.Results.solverName;
runID = parser.Results.runID;
outputFolder = parser.Results.outputFolder;
outputFileName = parser.Results.outputFileName;
printExcel = parser.Results.printExcel;
printText = parser.Results.printText;
printReport = parser.Results.printReport;
keepInputs = parser.Results.keepInputs;
keepGamsOutputs = parser.Results.keepGamsOutputs;
verbose = parser.Results.verbose;

% correct size of constrOpt
if ~isempty(constrOpt.rxnList)
    if size(constrOpt.rxnList, 1) > size(constrOpt.rxnList,2); constrOpt.rxnList = constrOpt.rxnList'; end;
    if size(constrOpt.values, 1) > size(constrOpt.values,2); constrOpt.values = constrOpt.values'; end;
end

% first, verify that GAMS is installed in your system
gamsPath = which('gams');
if isempty(gamsPath); error('OptForce: GAMS is not installed in your system. Please install GAMS.'); end;

%name of the function to solve the optimization problem in GAMS
gamsMustULFunction = 'findMustUL.gms';
%path of that function
pathGamsFunction = which(gamsMustULFunction);
if isempty(pathGamsFunction); error(['OptForce: ' gamsMustULFunction ' not in MATLAB path.']); end;
%current path
workingPath = pwd;
%go to the path associate to the ID for this run.
if ~isdir(runID); mkdir(runID); end; cd(runID);

% if the user wants to generate a report.
if printReport
    %create name for file.
    hour = clock;
    reportFileName = ['report-' date '-' num2str(hour(4)) 'h' '-' num2str(hour(5)) 'm.txt'];
    freport = fopen(reportFileName, 'w');
    reportClosed = 0;
    % print date of running.
    fprintf(freport, ['findMustULWithGAMS.m executed on ' date ' at ' num2str(hour(4)) ':' num2str(hour(5)) '\n\n']);
    % print matlab version.
    fprintf(freport, ['MATLAB: Release R' version('-release') '\n']);
    % print gams version.
    fprintf(freport, ['GAMS: ' regexprep(gamsPath, '\\', '\\\') '\n']);
    % print solver used in GAMS to solve optForce.
    fprintf(freport, ['GAMS solver: ' solverName '\n']);

    %print each of the inputs used in this running.
    fprintf(freport, '\nThe following inputs were used to run OptForce: \n');
    fprintf(freport, '\n------INPUTS------\n');
    %print model.
    fprintf(freport, '\nModel:\n');
    for i = 1:length(model.rxns)
        rxn = printRxnFormula(model, model.rxns{i}, false);
        fprintf(freport, [model.rxns{i} ': ' rxn{1} '\n']);
    end
    %print lower and upper bounds, minimum and maximum values for each of
    %the reactions in wild-type and mutant strain
    fprintf(freport, '\nLB\tUB\tMin_WT\tMax_WT\n');
    for i = 1:length(model.rxns)
        fprintf(freport, '%6.4f\t%6.4f\t%6.4f\t%6.4f\n', model.lb(i), model.ub(i), minFluxesW(i), maxFluxesW(i));
    end

    %print constraints
    fprintf(freport,'\nConstrained reactions:\n');
    for i = 1:length(constrOpt.rxnList)
        fprintf(freport,'%s: fixed in %6.4f\n', constrOpt.rxnList{i}, constrOpt.values(i));
    end

    fprintf(freport, '\nExcluded Reactions:\n');
    for i = 1:length(excludedRxns)
        rxn = printRxnFormula(model, excludedRxns{i}, false);
        fprintf(freport, [excludedRxns{i} ': ' rxn{1} '\n']);
    end

    fprintf(freport, '\nReactions from first order sets(MustU and MustL):\n');
    for i = 1:length(mustSetFirstOrder)
        rxn = printRxnFormula(model, mustSetFirstOrder{i}, false);
        fprintf(freport, [mustSetFirstOrder{i} ': ' rxn{1} '\n']);
    end

    fprintf(freport,'\nrunID(Main Folder): %s \n\noutputFolder: %s \n\noutputFileName: %s \n',...
        runID, outputFolder, outputFileName);

    fprintf(freport,'\nprintExcel: %1.0f \n\nprintText: %1.0f \n\nprintReport: %1.0f \n\nkeepInputs: %1.0f  \n\nkeepGamsOutputs: %1.0f \n\nverbose: %1.0f \n',...
        printExcel, printText, printReport, keepInputs, keepGamsOutputs, verbose);

end

copyfile(pathGamsFunction);

% export inputs for running the optimization problem in GAMS to find the
% MustUL Set
inputFolder = 'InputsMustUL';
exportInputsMustOrder2ToGAMS(model, 'UL', minFluxesW, maxFluxesW, constrOpt, excludedRxns, mustSetFirstOrder, inputFolder)

% create a directory to save results if this don't exist
if ~exist(outputFolder, 'dir')
    mkdir(outputFolder);
end

%run
if verbose
    run = system(['gams ' gamsMustULFunction ' lo=3 --myroot=' inputFolder '/ --solverName=' solverName ' gdx=GtoMUL']);
else
    run=system(['gams ' gamsMustULFunction ' --myroot=' inputFolder '/ --solverName=' solverName ' gdx=GtoMUL']);
end

if printReport; fprintf(freport, '\n------RESULTS------\n'); end;

%if user decide not to show inputs files for findMustUL.gms
if ~keepInputs; rmdir(inputFolder, 's'); end;

%if findMustUL.gms was executed correctly "run" should be 0
if run == 0

    if printReport; fprintf(freport, '\nGAMS was executed correctly\n'); end;
    if verbose; fprintf('GAMS was executed correctly\nSummary of information exported by GAMS:\n'); end;
    %show GAMS report in MATLAB console
    if verbose; gdxWhos GtoMUL; end;
    try
        findMustUL.name = 'findMustUL';
        rgdx('GtoMUL', findMustUL); %if do not exist the variable findMustUL in GtoMUL, an error will ocurr.
        if printReport; fprintf(freport, '\nGAMS variables were read by MATLAB correctly\n'); end;
        if verbose; fprintf('GAMS variables were read by MATLAB correctly\n'); end;

        %Using GDXMRW to read solutions found by findMustUL.gms
        %extract matrix 1 found by findMustII.gms. This matrix contains the
        %first reaction in each couple of reactions
        m1.name = 'matrix1';
        m1.compress = 'true';
        m1 = rgdx('GtoMUL', m1);
        uels_m1 = m1.uels{2};

        if ~isempty(uels_m1)
            %if the uel array for m1 is not empty, at least 1 couple of reations was found.
            if printReport; fprintf(freport, '\na MustUL set was found\n'); end;
            if verbose; fprintf('a MustUL set was found\n'); end;

            %find values for matrix 1
            val_m1 = m1.val;
            m1_full = full(sparse(val_m1(:,1), val_m1(:,2:end-1), val_m1(:,3)));

            %find values for matrix 2
            m2.name = 'matrix2';
            m2.compress = 'true';
            m2 = rgdx('GtoMUL', m2);
            uels_m2 = m2.uels{2};
            val_m2 = m2.val;
            m2_full = full(sparse(val_m2(:,1), val_m2(:,2:end-1), val_m2(:,3)));

            %initialize empty array for storing
            n_mustSet = size(m1_full,1);
            mustUL = cell(n_mustSet, 2);
            pos_mustUL = zeros(size(mustUL));
            mustUL_linear = {};

            %write each couple of reactions.
            for i = 1:n_mustSet
                rxn1 = uels_m1(m1_full(i,:) == 1);
                rxn2 = uels_m2(m2_full(i,:) == 1);
                mustUL(i,1) = rxn1;
                mustUL(i,2) = rxn2;
                pos_mustUL(i,1) = find(strcmp(model.rxns, rxn1));
                pos_mustUL(i,2) = find(strcmp(model.rxns, rxn2));
                mustUL_linear = union(mustUL_linear, [rxn1;rxn2]);
            end
            pos_mustUL_linear = cell2mat(arrayfun(@(x)find(strcmp(x, model.rxns)), mustUL_linear, 'UniformOutput', false))';
        else
            %if the uel array for m1 is empty, no couple of reations was found.
            if printReport; fprintf(freport, '\na MustUL set was not found\n'); end;
            if verbose; fprintf('a MustUL set was not found\n'); end;

            %initialize arrays to be returned by this function
            mustUL = {};
            pos_mustUL = [];
            mustUL_linear = {};
            pos_mustUL_linear = [];
        end

        % print info into an excel file if required by the user
        if printExcel
            if  ~isempty(uels_m1)
                currentFolder = pwd;
                cd(outputFolder);
                must = cell(size(mustUL,1), 1);
                for i = 1:size(mustUL, 1)
                    must{i} = strjoin(mustUL(i,:), ' or ');
                end
                xlswrite([outputFileName '_Info'],[{'Reactions'};must]);
                xlswrite(outputFileName, mustUL_linear);
                cd(currentFolder);
                if verbose
                    fprintf(['MustUL set was printed in ' outputFileName '.xls  \n']);
                    fprintf(['MustUL set was also printed in ' outputFileName '_Info.xls  \n']);
                end
                if printReport
                    fprintf(freport, ['\nMustUL set was printed in ' outputFileName '.xls  \n']);
                    fprintf(freport, ['\nMustUL set was printed in ' outputFileName '_Info.xls  \n']);
                end
            else
                if verbose; fprintf('No mustUL set was found. Therefore, no excel file was generated\n'); end;
                if printReport; fprintf(freport, '\nNo mustUL set was found. Therefore, no excel file was generated\n'); end;
            end
        end

        % print info into a plain text file if required by the user
        if printText
            if ~isempty(uels_m1)
                currentFolder = pwd;
                cd(outputFolder);
                f = fopen([outputFileName '_Info.txt'], 'w');
                fprintf(f,'Reactions\n');
                for i = 1:size(mustUL,1)
                    fprintf(f, '%s or %s\n', mustUL{i,1}, mustUL{i,2});
                end
                fclose(f);

                f = fopen([outputFileName '.txt'], 'w');
                for i = 1:length(mustUL_linear)
                    fprintf(f, '%s\n', mustUL_linear{i});
                end
                fclose(f);

                cd(currentFolder);
                if verbose
                    fprintf(['MustUL set was printed in ' outputFileName '.txt  \n']);
                    fprintf(['MustUL set was also printed in ' outputFileName '_Info.txt  \n']);
                end
                if printReport
                    fprintf(freport, ['\nMustUL set was printed in ' outputFileName '.txt  \n']);
                    fprintf(freport, ['\nMustUL set was printed in ' outputFileName '_Info.txt  \n']);
                end

            else
                if verbose; fprintf('No mustUL set was found. Therefore, no plain text file was generated\n'); end;
                if printReport; fprintf(freport, '\nNo mustUL set was found. Therefore, no plain text file was generated\n'); end;
            end
        end

        %close file for saving report
        if printReport; fclose(freport); reportClosed = 1; end;
        if printReport; movefile(reportFileName, outputFolder); end;
        delete(gamsMustULFunction);

        %remove or move additional files that were generated during running
        if keepGamsOutputs
            if ~isdir(outputFolder); mkdir(outputFolder); end;
            movefile('GtoMUL.gdx', outputFolder);
            movefile(regexprep(gamsMustULFunction, 'gms', 'lst'), outputFolder);
        else
            delete('GtoMUL.gdx');
            delete(regexprep(gamsMustULFunction, 'gms', 'lst'));
        end

        %go back to the original path
        cd(workingPath);
    catch
        %GAMS variables were not read correctly by MATLAB
        if verbose; fprintf('GAMS variables were not read by MATLAB corretly\n'); end;
        if printReport && ~reportClosed; fprintf(freport, '\nGAMS variables were not read by MATLAB corretly\n'); fclose(freport); end;
        cd(workingPath);
        error('OptForce: GAMS variables were not read by MATLAB corretly');

    end

    %if findMustUL.gms was not executed correctly "run" should be different from 0
else
    %if GAMS was not executed correcttly
    if printReport && ~reportClosed; fprintf(freport, '\nGAMS was not executed correctly\n'); fclose(freport); end;
    if verbose; fprintf('GAMS was not executed correctly\n'); end;
    cd(workingPath);
    error('OptForce: GAMS was not executed correctly');

end

end

author: opencobra

repo: cobratoolbox


text:
filename='Gripping_quad_fine';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'nonadjoint_compliance'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
constraint_case = 'INEQUALITY';
optimizer = 'MMA'; incrementFactor = 1;
filterType = 'P1';

nsteps = 1;
Vfrac_final = 0.5;
optimality_final = 1e-3;
constr_final =1e-3;

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

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function SO3F = cat(dim,varargin)
% overloads cat

% remove empty arguments
varargin(cellfun('isempty',varargin)) = [];
SO3F = varargin{1};
bw = SO3F.bandwidth;

for i = 2:numel(varargin)

  bw2 = varargin{i}.bandwidth;
  if bw > bw2
    varargin{i}.bandwidth = bw;
  else
    bw = bw2;
    SO3F.bandwidth = bw2;
  end

  ensureCompatibleSymmetries(SO3F,varargin{i});

  SO3F.fhat = cat(1+dim, SO3F.fhat, varargin{i}.fhat);
end

author: mtex-toolbox

repo: mtex

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


text:
function c8_le_li_test ( )

%*****************************************************************************80
%
%% C8_LE_LI_TEST tests C8_LE_LI.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'C8_LE_LI_TEST\n' );
  fprintf ( 1, '  C8_LE_LI evalues (C1 <= C2) using the Loo norm.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '        C1=C8_UNIFORM_01          C2=C8_UNIFORM_01         L3=C8_LE_LI(C1,C2)\n' );
  fprintf ( 1, '     ---------------------     ---------------------     ---------------------\n' );
  fprintf ( 1, '\n' );

  seed = 123456789;

  for test = 1 : 10

    [ c1, seed ] = c8_uniform_01 ( seed );
    [ c2, seed ] = c8_uniform_01 ( seed );
    l3 = c8_le_li ( c1, c2 );

    fprintf ( 1, '  %12.4f,%12.4f  %12.4f,%12.4f          %d\n', ...
      real ( c1 ), imag ( c1 ), real ( c2 ), imag ( c2 ), l3 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m