[Ricorsione] Divide & Conquer su un vettore

[elgiovo]

Ciao ragazzi,

stavo scrivendo un programmino Matlab (ma è un argomento trasversale che va bene per quasi tutti i linguaggi) che splitta ricorsivamente un vettore. La mia difficoltà sta nel fatto che il mio algoritmo non è una semplice ricerca binaria, che "butta via" una parte del vettore e tiene l'altra. Ho infatti una funzione piuttosto complicata che agisce su un vettore e trova al più un indice in cui si verifica una certa condizione. Per fare un esempio, la funzione può essere definita come

f(v) = indice(max(sqrt(v)))

che può trovare al più un indice (nessuno se v contiene solo numeri negativi).
Ora, se la funzione non trova nulla il codice deve fermarsi, mentre se trova un indice deve proseguire ad applicare la funzione ai due sottovettori che derivano dallo splitting del vettore originario. Lo splitting va effettuato nell'indice appena trovato. Il risultato finale dovrebbe essere un vettore con gli indici del vettore di partenza tali per cui la funzione f() ha dato un risultato. La parte in neretto mi sta creando delle difficoltà, perché ogni volta che vengono generati due sottovettori i loro indici partono sempre da 1. Avevo pensato di "incollare" sotto al vettore di partenza un vettore di indici e far poi agire la funzione sulla prima riga della matrice che ne risulta, ma non mi sembra molto elegante.

Siccome penso che sia un problema già affrontato da altri, vorre evitare di reinventare la ruota e vi chiedo umilmente aiuto.

Grazie!

[apatriarca]

Per implementare la somma cumulativa potevi anche fare uso della funzione cumsum di matlab. Per quanto riguarda le due righe del tipo:

if length(right) > 1
    rightIndices = step_detect(right,offset + i + 1);
else rightIndices = [];
end

Verificherei il caso all'interno della funzione e non prima di effettuare la ricorsione per maggiore chiarezza. In effetti vedo ora che dovrei aggiungere il test anche nel mio codice. Il codice mi sembra in effetti corretto comunque. Immagino tu non abbia provato a confrontare il tuo risultato con il mio codice.. Se riesco ci provo.

[elgiovo]

No non ho provato, è dai tempi di Informatica I (7 anni fa) che non uso C (lo so, sono deprecabile)
Se ti va prova, ma non pretendo che tu lo faccia, ci mancherebbe.

Nel frattempo ieri ho trovato anche questo toolbox di Matlab. Ho inviato richiesta per averlo, vediamo se mi rispondono.

[apatriarca]

Dopo qualche correzione del codice è riuscito a tirare fuori qualcosa:

  1. Change-point at index  1879 with confidence 100.00%
  2. Change-point at index   323 with confidence 100.00%
  3. Change-point at index   178 with confidence 100.00%
  4. Change-point at index    63 with confidence 100.00%
  5. Change-point at index    22 with confidence  99.93%
  6. Change-point at index   131 with confidence  98.46%
  7. Change-point at index   254 with confidence 100.00%
  8. Change-point at index  1272 with confidence 100.00%
  9. Change-point at index   957 with confidence 100.00%
 10. Change-point at index   725 with confidence 100.00%
 11. Change-point at index   534 with confidence 100.00%
 12. Change-point at index   407 with confidence 100.00%
 13. Change-point at index   344 with confidence  99.79%
 14. Change-point at index   331 with confidence  96.83%
 15. Change-point at index   597 with confidence 100.00%
 16. Change-point at index   678 with confidence  99.47%
 17. Change-point at index   705 with confidence  91.18%
 18. Change-point at index   851 with confidence 100.00%
 19. Change-point at index  1060 with confidence  99.08%
 20. Change-point at index   990 with confidence  95.87%
 21. Change-point at index   972 with confidence  98.89%
 22. Change-point at index  1015 with confidence  97.36%
 23. Change-point at index  1201 with confidence  91.21%
 24. Change-point at index  1196 with confidence  98.14%
 25. Change-point at index  1561 with confidence 100.00%
 26. Change-point at index  3134 with confidence 100.00%
 27. Change-point at index  2237 with confidence 100.00%
 28. Change-point at index  2171 with confidence 100.00%
 29. Change-point at index  1913 with confidence 100.00%
 30. Change-point at index  2115 with confidence  99.82%
 31. Change-point at index  2037 with confidence  99.96%
 32. Change-point at index  2197 with confidence  98.60%
 33. Change-point at index  2176 with confidence  99.32%
 34. Change-point at index  2787 with confidence 100.00%
 35. Change-point at index  2493 with confidence 100.00%
 36. Change-point at index  2287 with confidence 100.00%
 37. Change-point at index  3347 with confidence 100.00%
 38. Change-point at index  3240 with confidence  99.46%
 39. Change-point at index  3202 with confidence  99.94%
 40. Change-point at index  3155 with confidence  99.30%
 41. Change-point at index  3317 with confidence  99.84%
 42. Change-point at index  3537 with confidence 100.00%
 43. Change-point at index  3430 with confidence 100.00%
 44. Change-point at index  3383 with confidence  96.10%
 45. Change-point at index  3355 with confidence  90.46%
 46. Change-point at index  3418 with confidence  93.63%
 47. Change-point at index  3636 with confidence  98.74%

Il nuovo codice:

#include <float.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>


/* Recursive change point analysis using CUSUM estimator */
unsigned change_point_analysis(double *data_in, unsigned size, unsigned offset, unsigned *indices_out, double *confidences)
{
    if (size == 0) { return 0; }

    /* Finds average */
    double average = 0.0;
    for (unsigned i = 0; i < size; ++i) {
        average += data_in[i];
    }
    average /= size;

    /* Finds minimum and maximum of CUSUM chart */
    double S = 0.0, Smin = 0.0, Smax = 0.0, Sdiff = 0.0;
    int max_index = -1;
    int min_index = -1;
    for (unsigned i = 0; i < size; ++i) {
        S += data_in[i] - average;
        if (S < Smin) {
            Smin = S;
            min_index = (int)i;
        }
        if (S > Smax) {
            Smax = S;
            max_index = (int)i;
        }
    }
    Sdiff = Smax - Smin;
    unsigned index = fabs(Smin) < fabs(Smax) ? max_index : min_index;

    /* Bootstrapping */
    srand((unsigned int)time(NULL));

    double * bootstrap_array = (double *)malloc(sizeof(double)*size);
    if (!bootstrap_array) {
        fputs("Unable to allocate memory.\n", stderr);
        exit(1);
    }

    memcpy(bootstrap_array, data_in, sizeof(double)*size);

    double confidence = 0.0;
    for (int k = 0; k < 10000; ++k) {
        /* Random permutation */
        for (unsigned i = 0; i < size; ++i) {
            int j = rand() % (i + 1);
            double tmp = bootstrap_array[i];
            bootstrap_array[i] = bootstrap_array[j];
            bootstrap_array[j] = tmp;
        }

        /* CUSUM chart computation */
        double Sb = 0.0, Sbmin = 0.0, Sbmax = 0.0, Sbdiff = 0.0;
        for (unsigned i = 0; i < size; ++i) {
            Sb += bootstrap_array[i] - average;
            Sbmin = Sb < Sbmin ? Sb : Sbmin;
            Sbmax = Sb > Sbmax ? Sb : Sbmax;
        }
        Sbdiff = Sbmax - Sbmin;
        if (Sbdiff < Sdiff) { ++confidence; }
    }
    confidence /= 10000.0;

    free(bootstrap_array);
    if (confidence < 0.9) { /* Voglio confidenza di almeno il 90% */
        return 0;
    }

    /* Output and recursion */
    indices_out[0] = index + offset;
    confidences[0] = confidence;
    unsigned ret = 1;
    ret += change_point_analysis(data_in, index, offset, indices_out+ret, confidences+ret);
    ret += change_point_analysis(data_in+index+1, size-index-1, offset+index+1, indices_out+ret, confidences+ret);
    return ret;
}

int main(void)
{
    double *data = NULL;
    unsigned size = 0;
    unsigned capacity = 0;
    unsigned *change_points = NULL;
    double *confidences = NULL;

    FILE * file = fopen("timeseries.txt", "r");
    if (!file) {
        fputs("Unable to open input file.\n", stderr);
        return 1;
    }

    int i = -1;
    do {
        ++i;
        if (capacity <= i) {
            unsigned new_capacity = capacity > 0 ? 2*capacity : 16;
            double *new_data = (double *)malloc(sizeof(double)*new_capacity);
            if (!new_data) {
                fputs("Unable to allocate memory.\n", stderr);
                free(data);
                fclose(file);
                exit(1);
            }
            memcpy(new_data, data, i*sizeof(double));
            free(data);
            data = new_data;
            capacity = new_capacity;
        }
    } while (fscanf(file, "%f", data + i) > 0);
    size = i;

    fclose(file);

    change_points = (unsigned *)malloc(sizeof(unsigned)*size);
    if (!change_points) {
        fputs("Unable to allocate memory.\n", stderr);
        free(data);
        exit(1);
    }
    memset(change_points, 0, sizeof(unsigned)*size);

    confidences = (double *)malloc(sizeof(double)*size);
    if (!confidences) {
        fputs("Unable to allocate memory.\n", stderr);
        free(change_points);
        free(data);
        exit(1);
    }
    memset(confidences, 0, sizeof(double)*size);

    unsigned n = change_point_analysis(data, size, 1, change_points, confidences);
    for (unsigned i = 0; i < n; ++i) {
        printf("%3d. Change-point at index %5d with confidence %6.2f%%\n", i+1, change_points[i], confidences[i]*100.0);
    }

    free(data); free(change_points); free(confidences);
}

[elgiovo]

Ok! Questi sono gli indici che ha trovato il tuo programma.



Il mio programma trova il salto più grosso (1879) poi si impasta nel sottovettore sinistro dove comincia a trovare una sfilza di 1. Devo cercare di capire dov'è la differenza rispetto al tuo.

C'è poi da fare una taratura dei parametri, tipo il numero di bootstrap e la soglia per l'intervallo di confidenza, perché i salti rilevati sono troppi. Ad esempio, tra gli indici 2000 e 4000 io metterei non più di due salti (i più evidenti).

Intanto grazie per l'interesse dimostrato!

[apatriarca]

Se vuoi ridurre il numero di punti, forse conviene anche aggiungere un massimo livello di ricorsione. Anche aumentando il livello di confidenza sopra il 95% sarebbero infatti ancora molti i valori estratti. Ma l'articolo mi sembra parlasse anche di metodi più avanzati per l'eliminazione dei falsi positivi.

[elgiovo]

"apatriarca":Se vuoi ridurre il numero di punti, forse conviene anche aggiungere un massimo livello di ricorsione. Anche aumentando il livello di confidenza sopra il 95% sarebbero infatti ancora molti i valori estratti. Ma l'articolo mi sembra parlasse anche di metodi più avanzati per l'eliminazione dei falsi positivi.

Ciao scusa ma non ho più avuto tempo da dedicare al codice.
In effetti mi pare che aumentare o diminuire quel 95% non sortisca effetti... Mah..
Il problema che avevo con gli "1" trovati dal codice era nello stimatore dell'indice del salto. Io usavo quello con la varianza, mentre il problema va via se uso quello che deriva dalle somme cumulative (mi pare che sia quello che usi nel tuo impeccabile codice C).
Si mi ricordo che parlava dell'eliminazione di falsi positivi, se trovo tempo approfondisco.
Il metodo comunque è piuttosto interessante (come vedi porta a pubblicazioni), però ho paura che i miei dati siano un pò troppo rumorosi.
Limitando il massimo della ricorsione i risultati non sono male, però continua a trovarmi salti dove a occhio non li metterei e a non metterli dove a occhio li metterei. A tempo perso ci lavoro, se ho aggiornamenti li condivido.

[apatriarca]

So che probabilmente è un'idea un po' naïve ma motivata da alcuni degli esami fatti ad ingegneria.. Perché non filtri il segnale con un filtro passa basso che elimina (attenua) le componenti di alta frequenza?