[Vari] Contest 8: Compattamento Barionico

[gugoXX]

1) Sia data una lista di punti contenuti in uno spazio cubico, con coordinate tridimensionali in virgola mobile.
Il formato della lista e'

NP Size
X1 Y1 Z1
X2 Y2 Z2
...
Xn Yn Zn

Dove NP e' il numero di punti (quindi le linee che seguono)
e Size e' un numero floating point pari alla dimensione del cubo di contenimento (quindi nessuna coordinata superera' tale valore)

Trovare ed enunciare la coppia di punti a distanza minima e la coppia di punti a distanza massima (e le distanze relative)

PS: Per quanto riguarda la formula della distanza tridimensionale in uno spazio euclideo, dati i punti
Pa di coordinate (Xa, Ya, Za) e Pb di coordinate (Xb, Yb, Zb)
http://operaez.net/mimetex/dist(Pa,Pb)=sqrt((Xa-Xb)%5E2+(Ya-Yb)%5E2+(Za-Zb)%5E2)

(Da Wiki, sotto distanza Euclidea: http://upload.wikimedia.org/math/2/e/2/2e2f050e0065e30cac6f0bd5df46a1fd.png )
E non ditemi che vi arenate sulla distanza 3D per favore. E' pieno di gente qui che vorrebbe fare giochi 3D, e chi si arena su una semplice distanza...

Lista dei punti: http://www.usaupload.net/d/cjepq9exi1j

2) Teorema: Sia dato un cubo di lato d ed un numero di particelle N
Esiste almeno una configurazione spaziale casuale delle N particelle tale per cui non esista nemmeno una coppia di particelle
con distanza c inferiore a
d / ³√N


http://img.photobucket.com/albums/v314/gugogugo/mimetex.gif

La semplice dimostrazione è lasciata come esercizio al lettore (Ho sempre sognato di scriverlo)

In pratica questo teorema ci assicura che tirando a caso molte volte la posizione di ciascuna delle N particelle, c'e' almeno una configurazione tale per cui le particelle non sono troppo vicine e distano tra loro almeno c.

Problema:
Si trovi almeno una configurazione spaziale (casuale) delle particelle nei seguenti casi:
d = 10 N = 5
d = 10 N = 10
d = 10 N = 50
d = 10 N = 100
d = 10 N = 200

PS: La soluzione di questo problema è in uso al CERN, per le condizioni iniziali di alcune simulazioni.

Domani posterò anche l'algoritmo di confronto a forza bruta.
Il primo esercizio serve come test di verifica per il secondo, non mi aspetto un'ottimizzazione particolare, anche se magari potrebbe essere interessante.

[banryu79]

Il primo esercizio serve come test di verifica per il secondo, non mi aspetto un'ottimizzazione particolare, anche se magari potrebbe essere interessante.
Dato che quello che interessa non è la reale distanza 3D ma solo il confronto di valori riguardanti le varie distanze dei vari punti per trovare la coppia con la minima distanza e la coppia con la massima distanza, nella formula per calcolare la distanza di un punto non dovrebbe essere neccessario fare la radice quadrata, si può così evitarne la sua chiamata.
(Sì lo so che parlavi di ottimizzazioni a livello algoritmico, però mi è saltata all'occhio questa cosa).

[gugoXX]

Dai cos'e'? Passiamo dalle cose quadrate a quelle sferiche e ci "Impalliamo" subito?

[banryu79]

Mah, non saprei, non è che sia perchè il Forum è un po' smorto in quest'ultimo periodo?
Magari gioverebbe anche capire meglio che cosa si sia alla fine deciso/accolto nella discussione circa i nuovi contest... link (http://www.hwupgrade.it/forum/showthread.php?t=1872715)

[Tommo]

Non capisco... ma per "casuale" si intende che va bene una configurazione qualsiasi, oppure che i valori iniziali devono essere tirati a caso?
C'è una bella differenza...

Nel primo caso la soluzione più veloce è:
nextPos = pos + c;
Sarebbe ben facile :asd:

Nel secondo invece c'è un problema... la configurazione casuale migliore potrebbe uscire alla seconda iterazione come dopo 1000... rendendo difatti insensato il contest...
credo sia il caso di fare un file che contenga le varie configurazioni iniziali.

[Vincenzo1968]

Ohé,

per me che, in matematica, sono deboluccio: è corretto il seguente output?:

http://www.guidealgoritmi.it/images/ImgForums/Contest08Soluz01.jpg

Questo è il codice con l'array di prova:


#include <stdio.h>
#include <math.h>

typedef struct tagPunti
{
double A;
double B;
double C;
} Punti;

typedef struct tagRes
{
Punti P1;
Punti P2;
double Distanza;
int index1;
int index2;
} Res;

void Calcola(Punti P[], int dim, Res *pMin, Res *pMax)
{
int x, y;
Punti P1, P2;
double distanza, distanzaMin, distanzaMax;

distanzaMin = distanzaMax = -1;

x = y = 0;

while( x < dim - 1 )
{
P1 = P[x];
y = x + 1;
while( y < dim )
{
P2 = P[y];
distanza = sqrt( pow(P1.A - P2.A, 2.0) + pow(P1.B - P2.B, 2.0) + pow(P1.C - P2.C, 2.0) );
//printf("Distanza -> %.15lf\n", distanza);
if ( distanza < distanzaMin || distanzaMin == -1 )
{
distanzaMin = distanza;
pMin->P1 = P[x];
pMin->P2 = P[y];
pMin->Distanza = distanza;
pMin->index1 = x;
pMin->index2 = y;
}
if ( distanza > distanzaMax )
{
distanzaMax = distanza;
pMax->P1 = P[x];
pMax->P2 = P[y];
pMax->Distanza = distanza;
pMax->index1 = x;
pMax->index2 = y;
}
++y;
}
++x;
}
}

int main()
{
Punti P[5];
Res r1, r2;

P[0].A = 1.5;
P[0].B = 2.3;
P[0].C = 5.1;

P[1].A = 2.8;
P[1].B = 4.3;
P[1].C = 7.5;

P[2].A = 8.5;
P[2].B = 2.3;
P[2].C = 4.1;

P[3].A = 5.5;
P[3].B = 8.1;
P[3].C = 6.1;

P[4].A = 7.5;
P[4].B = 9.3;
P[4].C = 2.6;

Calcola(P, 5, &r1, &r2);

printf("\nDistanza Min: P[%d] - P[%d] -> %.15lf\n", r1.index1, r1.index2, r1.Distanza);
printf("\nDistanza Max: P[%d] - P[%d] -> %.15lf\n", r2.index1, r2.index2, r2.Distanza);

return 0;
}

[gugoXX]

Non capisco... ma per "casuale" si intende che va bene una configurazione qualsiasi, oppure che i valori iniziali devono essere tirati a caso?
C'è una bella differenza...

Nel primo caso la soluzione più veloce è:
nextPos = pos + c;
Sarebbe ben facile :asd:

Nel secondo invece c'è un problema... la configurazione casuale migliore potrebbe uscire alla seconda iterazione come dopo 1000... rendendo difatti insensato il contest...
credo sia il caso di fare un file che contenga le varie configurazioni iniziali.

Casuale perche' o riesci a trovare una soluzione a priori "cristallina" a questo problema, e non e' semplicissimo quando non si hanno numeri di particelle "belli", oppure mi sa tanto che quella da te proposta non e' sufficientemente esaustiva.
Mi spiego, a forza di aggiungere +c alla coordinata da te scelta, prima o poi cozzerai contro la scatola, e allora cosa farai?

L'idea sarebbe quella di trovare un metodo per la ricerca delle soluzioni casuali, non un'unica soluzione deterministica e sempre uguale a fronte dei valori scelti, ma anche cosi' va bene, c'e' da divertirsi.

[Tommo]

Casuale perche' o riesci a trovare una soluzione a priori "cristallina" a questo problema, e non e' semplicissimo quando non si hanno numeri di particelle "belli", oppure mi sa tanto che quella da te proposta non e' sufficientemente esaustiva.
Mi spiego, a forza di aggiungere +c alla coordinata da te scelta, prima o poi cozzerai contro la scatola, e allora cosa farai?

L'idea sarebbe quella di trovare un metodo per la ricerca delle soluzioni casuali, non un'unica soluzione deterministica e sempre uguale a fronte dei valori scelti, ma anche cosi' va bene, c'e' da divertirsi.

vabè quella non era una soluzione seria :asd:

Era solo per dire che il problema non è affatto difficile se posso mettermi le particelle come mi pare... basta pensarci un pò di più per capire che per occupare il volume minore devi affiancare tutti tetraedroni, come nei cristalli... in quella maniera ogni particella è al centro di un poliedro di raggio c... in 2d sarebbero tutti esagoni.

E' un pò difficile a farsi ma credo non si possano compattare di più...

[marco.r]

Casuale perche' o riesci a trovare una soluzione a priori "cristallina" a questo problema, e non e' semplicissimo quando non si hanno numeri di particelle "belli", oppure mi sa tanto che quella da te proposta non e' sufficientemente esaustiva.
Mi spiego, a forza di aggiungere +c alla coordinata da te scelta, prima o poi cozzerai contro la scatola, e allora cosa farai?

Continuo con la fila successiva :D. E da quattro conti che ho fatto sembra funzionare per le varie scelte di N.

In un cubo mi ci stanno ( floor(d/c)+1 )^3 sfere (nell'ipotesi che una molecola possa stare "appicicata al bordo").
Quindi diciamo almeno (d/c)^3
Sostituendo la definizione di c ottieni (d/d * pow(N,1/3))^3 molecole, che guardacaso coincide proprio con N.
Dimostra contemporaneamente il teorema sopra e mostra un modo per impacchettare le molecole (anzi, scommetto che e' il modo in cui e' stato ricavato quel valore).
Quindi la soluzione "cristallina" e' effettivamente banale.

[gugoXX]

Esatto! Questa e' la dimostrazione che pero' funziona per "cubi", nel senso quando le particelle sono 1,8,27, etc.

Quando sono... 17, impaccarle in quel modo non e' semplice (e mi sa che non funziona). Pero' il teorema e' sempre valido.

[Vincenzo1968]

Mbe'? Nessuno sa dirmi se ho eseguito correttamente i calcoli?

Data la seguente lista di punti:


P[0] -> 1.5 2.3 5.1
P[1] -> 2.8 4.3 7.5
P[2] -> 8.5 2.3 4.1
P[3] -> 5.5 8.1 6.1
P[4] -> 7.5 9.3 2.6


è corretto dire che la coppia di punti a distanza minima è (P[0],P[1]) e quella a distanza massima è, invece, (P[0],P[4])?

[marco.r]

vabè quella non era una soluzione seria :asd:

Era solo per dire che il problema non è affatto difficile se posso mettermi le particelle come mi pare... basta pensarci un pò di più per capire che per occupare il volume minore devi affiancare tutti tetraedroni, come nei cristalli... in quella maniera ogni particella è al centro di un poliedro di raggio c... in 2d sarebbero tutti esagoni.

E' un pò difficile a farsi ma credo non si possano compattare di più...

Non penso sia quella ottimale, resta molto spazio tra le sfere.
Si potrebbe provare a vedere quante sfere si rescono effettivamente a compattare al massimo

[Tommo]

Non penso sia quella ottimale, resta molto spazio tra le sfere.
Si potrebbe provare a vedere quante sfere si rescono effettivamente a compattare al massimo

E' ottimale :D

Semplicemente, TUTTE le coppie di sfere sono distanti c, che è il minimo.
E se tutte quante le sfere sono fra loro distanti c, non si possono compattare di più...

prendi il Geomag per esempio: la struttura più densa possibile è quella composta da piramidi a base triangolare... se uno prende le palline come particelle, e le stecche come la distanza minima c.
Non c'è maniera di compattarlo di più, e credo valga lo stesso per i barioni, qualunque cosa siano :asd:

[Torav]

E' ottimale :D

Semplicemente, TUTTE le coppie di sfere sono distanti c, che è il minimo.
E se tutte quante le sfere sono fra loro distanti c, non si possono compattare di più...

prendi il Geomag per esempio: la struttura più densa possibile è quella composta da piramidi a base triangolare... se uno prende le palline come particelle, e le stecche come la distanza minima c.
Non c'è maniera di compattarlo di più, e credo valga lo stesso per i barioni, qualunque cosa siano :asd:

In verità i reticoli con il packing più alto (74.0% del volume occupato dalle pallozze :p ) sono diversi: ad esempio quello a base esagonale e il reticolo cubico a facce centrate. Per convincersi della cosa basta pensare di "stendere" delle palline su un piano in modo che siano il più vicino possibile. Ora mettiamo il secondo strato, le palline finiranno proprio negli avvallamenti del primo strato. Se ora mettiamo un terzo strato nella stessa posizione del primo otteniamo un reticolo esagonale, se invece lo mettiamo ancora girato di 60 gradi otteniamo l'FCC. Se provate con le biglie si vede chiaramente che più strette di così non ci possono stare :p
se domani ho tempo voglio vedere cosa riesco a tirare fuori :)

[Tommo]

In verità i reticoli con il packing più alto (74.0% del volume occupato dalle pallozze :p ) sono diversi: ad esempio quello a base esagonale e il reticolo cubico a facce centrate. Per convincersi della cosa basta pensare di "stendere" delle palline su un piano in modo che siano il più vicino possibile. Ora mettiamo il secondo strato, le palline finiranno proprio negli avvallamenti del primo strato. Se ora mettiamo un terzo strato nella stessa posizione del primo otteniamo un reticolo esagonale, se invece lo mettiamo ancora girato di 60 gradi otteniamo l'FCC. Se provate con le biglie si vede chiaramente che più strette di così non ci possono stare :p
se domani ho tempo voglio vedere cosa riesco a tirare fuori :)

Molto interessante... ma mi sfugge una cosa: se ruoti un esagono di 60 gradi, non cambia niente...
Per cui non capisco come fa la rotazione a compattarlo di più :mbe:

[gugoXX]

Per stimolare i giovani virgulti ho aggiornato il punto1, migliorandone la definizione e estendendolo con un file di prova.
Ora cerco di mettere il primo risultato forza bruta di base.

[gugoXX]

Ecco

Max Dist: 1.70445415447397 Indexes: 48483-72832
Min Dist: 0.000366913641159708 Indexes: 22503-72620

234827 (Millisecondi)


Qui il codice.
Vi chiederei quanto piu' possibile di tenere separate la fase di input da quella di elaborazione, qualunque essa sia.
Il risultato della fase di input vorrebbe essere una semplice lista o array in memoria, da elaborare successivamente.

PS: C'e' gia' una prima banalissima ottimizzazione, come gia' proposto, ovvero di non confrontare le vere e proprie distanze, ma il loro quadrato
in modo tale da evitare tante, tantissime radici quadrate.
Nel file di prova si risparmiano la bellezza di 4999950000 radici quadrate.


class Program
{
static void Main(string[] args)
{
DistCalc dc = new DistCalc();
dc.Load(@"C:\temp\Coords.dat");
Stopwatch sw = new Stopwatch();
sw.Start();
var res = dc.Research();
sw.Stop();
Console.WriteLine(res);
Console.WriteLine(sw.ElapsedMilliseconds);
Console.ReadKey();
}
}

public class DistCalc
{
List<Point3d> allpoints;
public DistCalc()
{
allpoints = new List<Point3d>();
}

public void Load(string fname)
{
StreamReader sr = new StreamReader(fname);
string stt = sr.ReadLine();
string[] first = stt.Split(new char[] { ' ' }, StringSplitOptions.RemoveEmptyEntries);
int howmany = int.Parse(first[0]);
InputSize = double.Parse(first[1]);
string line;
while ((line = sr.ReadLine()) != null)
{
string[] coords = line.Split(new char[] { ' ' }, StringSplitOptions.RemoveEmptyEntries);
allpoints.Add(new Point3d(double.Parse(coords[0]),double.Parse(coords[1]),double.Parse(coords[2])));
}
sr.Close();
}

public void Save(string fname)
{
StreamWriter sw = new StreamWriter(fname);
sw.WriteLine("{0} {1}", allpoints.Count, InputSize);
foreach (Point3d p3d in allpoints)
{
sw.WriteLine(p3d.ToString());
}
sw.Close();
}

public double InputSize;
public void Generate(int npoints,double dim)
{
InputSize = dim;
allpoints.Clear();
for (int t = 0; t < npoints; t++)
{
allpoints.Add(Point3d.Rand(dim));
}
}

public class Result
{
public Result()
{
curdimmin = double.MaxValue;
curdimmax = 0d;
}

public double curdimmin;
public int i1dimmin, i2dimmin;

public double curdimmax;
public int i1dimmax, i2dimmax;

public override string ToString()
{
StringWriter sw = new StringWriter();
sw.WriteLine("Max Dist: {0} Indexes: {1}-{2}", Math.Sqrt(curdimmax), i1dimmax, i2dimmax);
sw.WriteLine("Min Dist: {0} Indexes: {1}-{2}", Math.Sqrt(curdimmin), i1dimmin, i2dimmin);
return sw.ToString();
}
}

public Result Research()
{
Result res = new Result();
int ln = allpoints.Count;
for (int t = 0; t < ln - 1; t++)
{
for (int u = t+1; u < ln; u++)
{
double cd = allpoints[t] - allpoints[u];
if (cd < res.curdimmin)
{
res.curdimmin = cd;
res.i1dimmin = t;
res.i2dimmin = u;
}
else if (cd > res.curdimmax)
{
res.curdimmax = cd;
res.i1dimmax = t;
res.i2dimmax = u;
}
}
}
return res;
}
}

public struct Point3d
{
public double x;
public double y;
public double z;
public Point3d(double ix,double iy,double iz)
{
x = ix;
y = iy;
z = iz;
}

public static Random rnd = new Random(18);
public static Point3d Rand(double dim)
{
return new Point3d(rnd.NextDouble() * dim, rnd.NextDouble() * dim, rnd.NextDouble() * dim);
}

public static double operator -(Point3d d1, Point3d d2)
{
double dx=d2.x-d1.x;
double dy=d2.y-d1.y;
double dz=d2.z-d1.z;
return (dx * dx + dy * dy + dz * dz);
}

public override string ToString()
{
return string.Format("{0} {1} {2}", x, y, z);
}
}

[Vincenzo1968]

Per stimolare i giovani virgulti ho aggiornato il punto1, migliorandone la definizione e estendendolo con un file di prova.
Ora cerco di mettere il primo risultato forza bruta di base.

Cosa rappresenta quell'uno dopo il diecimila nella prima riga del file? :confused:

[Torav]

Molto interessante... ma mi sfugge una cosa: se ruoti un esagono di 60 gradi, non cambia niente...
Per cui non capisco come fa la rotazione a compattarlo di più :mbe:

ups! Era tardi..quello che intendo dire è che il terzo strato si può mettere in due modi diversi riempiendo esattamente la stessa % di spazio! Ora sto uscendo, se trovo qualche disegno più chiaro lo posto perchè altrimenti ci vuole una bella dose di immaginazione!

[Torav]

ecco un link (http://www.tf.uni-kiel.de/matwis/amat/def_en/kap_5/backbone/r5_4_1.html) che fa vedere decentemente la cosa, e con questo chiudo il leggero OT, scusate :)

[Vincenzo1968]

La mia soluzione brute-force per il punto 1:

http://www.guidealgoritmi.it/images/ImgForums/Contest08ARis01.jpg

La mia macchina:

AMD Athlon(tm) 64 X2
Dual Core Processor 4800+
2.50 GHz
896 MB di RAM

Microsoft Windows XP Professional (32 bit)
Service Pack 3


El codigo fuente:

#include <stdio.h>
#include <math.h>
#include <time.h>
#include <malloc.h>

#define BUFFER_SIZE 4096

typedef struct tagPunti
{
double A;
double B;
double C;
} Punti;

typedef struct tagRes
{
Punti P1;
Punti P2;
double Distanza;
int index1;
int index2;
} Res;

typedef struct tagArray
{
int dimensione;
double dmax;
Punti *m;
} Array;

typedef enum tagStati
{
S_ERROR = -1, S0 = 0, S1, S2, S3, S4, S5, S6, S7, S8, S9
} Stati;


Stati DFA(const char *szFileName, Array *pArray)
{
Stati stato = S0;
FILE *fp;
unsigned char buffer[BUFFER_SIZE];
int numblocks;
int numread;
unsigned char c;
int k, x, j;
int riga, colonna;
char szNum[256];
double num;

unsigned char byteCR = 0xD; // Carriage Return
unsigned char byteLF = 0xA; // Line Feed

fp = fopen(szFileName, "rb");
if ( fp == NULL )
{
printf("Errore nell'apertura del file %s\n", szFileName);
return S_ERROR;
}

if ( fseek(fp, 0, SEEK_END) )
return 0;

numblocks = ftell(fp)/BUFFER_SIZE;
if ( numblocks == 0 )
{
numblocks = 1;
}
else
{
if ( ftell(fp) % BUFFER_SIZE != 0 )
numblocks++;
}

fseek(fp, 0, SEEK_SET);
numread = fread(buffer, 1, BUFFER_SIZE, fp);
if (numread == 0)
{
fclose(fp);
printf("Errore 1 nella lettura del file %s\n", szFileName);
return S_ERROR;
}

pArray->dimensione = 0;
k = 0;
c = *(buffer + k++);
while ( 1 )
{
if ( c >= '0' && c <= '9' )
{
pArray->dimensione = pArray->dimensione * 10 + c - '0';
}
else if ( c == ' ' || c == '\t' )
{
j = 0;
while ( c != byteLF )
{
c = *(buffer + k++);
if ( c >= '0' && c <= '9' )
szNum[j++] = c;
}
szNum[j] = '\0';
pArray->dmax = atof(szNum);
break;
}
else if ( c == '\n' )
{
break;
}

c = *(buffer + k++);
}

pArray->m = (Punti*)malloc(pArray->dimensione*sizeof(Punti));
if ( !pArray->m )
{
printf("Memoria insufficiente.\n");
fclose(fp);
return S_ERROR;
}

riga = colonna = 0;
num = 0;
x = j = 0;
while ( x < numblocks )
{
c = *(buffer + k++);

if ( c == byteLF )
{
riga++;
colonna = 0;
}

switch (stato)
{
case S0:
j = 0;
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S2;
}
else if ( c == '.' )
{
szNum[j++] = c;
stato = S3;
}
else if ( c == '-' || c == '+' )
{
szNum[j++] = c;
stato = S1;
}
else if ( c == ' ' || c == '\t' )
{
while ( c == ' ' || c == '\t' )
{
c = *(buffer + k++);
if ( k >= numread )
break;
}
k--;
}
else if ( c == byteCR || c == byteLF )
{
// nada
}
else
{
printf("Errore: stato S0; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S1:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S2;
}
else if ( c == '.' )
{
szNum[j++] = c;
stato = S3;
}
else
{
printf("Errore: stato S1; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S2:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
}
else if ( c == '.' )
{
szNum[j++] = c;
stato = S4;
}
else if ( c == 'E' || c == 'e' )
{
szNum[j++] = c;
stato = S5;
}
else if ( c == ' ' || c == '\t' || c == byteCR || c == byteLF )
{
szNum[j] = '\0';
if ( colonna == 0 )
pArray->m[riga].A = atof(szNum);
else if ( colonna == 1 )
pArray->m[riga].B = atof(szNum);
else if ( colonna == 2 )
pArray->m[riga].C = atof(szNum);
colonna++;
stato = S0;
}
else
{
printf("Errore: stato S2; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S3:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S6;
}
else
{
printf("Errore: stato S3; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S4:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S6;
}
else if ( c == 'E' || c == 'e' )
{
szNum[j++] = c;
stato = S7;
}
else if ( c == ' ' || c == '\t' || c == byteCR || c == byteLF )
{
szNum[j] = '\0';
if ( colonna == 0 )
pArray->m[riga].A = atof(szNum);
else if ( colonna == 1 )
pArray->m[riga].B = atof(szNum);
else if ( colonna == 2 )
pArray->m[riga].C = atof(szNum);
colonna++;
stato = S0;
}
else
{
printf("Errore: stato S4; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S5:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S9;
}
else if ( c == '-' || c == '+' )
{
szNum[j++] = c;
stato = S8;
}
else
{
printf("Errore: stato S5; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S6:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S6;
}
else if ( c == 'E' || c == 'e' )
{
szNum[j++] = c;
stato = S7;
}
else if ( c == ' ' || c == '\t' || c == byteCR || c == byteLF )
{
szNum[j] = '\0';
if ( colonna == 0 )
pArray->m[riga].A = atof(szNum);
else if ( colonna == 1 )
pArray->m[riga].B = atof(szNum);
else if ( colonna == 2 )
pArray->m[riga].C = atof(szNum);
colonna++;
stato = S0;
}
else
{
printf("Errore: stato S6; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S7:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S9;
}
else if ( c == '-' || c == '+' )
{
szNum[j++] = c;
stato = S8;
}
else
{
printf("Errore: stato S7; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S8:
case S9:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S9;
}
else if ( c == ' ' || c == '\t' || c == byteCR || c == byteLF )
{
szNum[j] = '\0';
if ( colonna == 0 )
pArray->m[riga].A = atof(szNum);
else if ( colonna == 1 )
pArray->m[riga].B = atof(szNum);
else if ( colonna == 2 )
pArray->m[riga].C = atof(szNum);
colonna++;
stato = S0;
}
else
{
printf("Errore: stato S9; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;
}

if ( k >= numread )
{
numread = fread(buffer, 1, BUFFER_SIZE, fp);
k = 0;
x++;
}
}

fclose(fp);

return stato;
}

void Calcola(Punti P[], int dim, Res *pMin, Res *pMax)
{
int x, y;
double distanza, distanzaMin, distanzaMax;
double diffA, diffB, diffC;

distanzaMin = distanzaMax = -1;

x = y = 0;

while( x < dim - 1 )
{
y = x + 1;
while( y < dim )
{
//distanza = sqrt( pow(P[x].A - P[y].A, 2.0) + pow(P[x].B - P[y].B, 2.0) + pow(P[x].C - P[y].C, 2.0) );
diffA = P[x].A - P[y].A;
diffB = P[x].B - P[y].B;
diffC = P[x].C - P[y].C;
distanza = diffA*diffA + diffB*diffB + diffC*diffC;

if ( distanza < distanzaMin || distanzaMin == -1 )
{
distanzaMin = distanza;
pMin->P1 = P[x];
pMin->P2 = P[y];
pMin->index1 = x;
pMin->index2 = y;
}
if ( distanza > distanzaMax )
{
distanzaMax = distanza;
pMax->P1 = P[x];
pMax->P2 = P[y];
pMax->index1 = x;
pMax->index2 = y;
}
++y;
}
++x;
}

pMin->Distanza = sqrt( distanzaMin );
pMax->Distanza = sqrt( distanzaMax );
}

int main()
{
Stati stato;
Array myArray;
Res r1, r2;

clock_t c_start, c_end;

char *szFileName = "C:\\Scaricamenti\\Temp\\Gugo\\Contest 08 - Compattamento Barionico\\Coords\\Coords.dat";

c_start = clock();
stato = DFA(szFileName, &myArray);
c_end = clock();
printf("\nTempo impiegato per caricamento dati -> %5.5f secondi\n", (double)(c_end - c_start) / CLOCKS_PER_SEC);

c_start = clock();
Calcola(myArray.m, myArray.dimensione, &r1, &r2);
c_end = clock();

printf("\nDistanza Min: P[%d] - P[%d] -> %.15lf\n", r1.index1, r1.index2, r1.Distanza);
printf("\nDistanza Max: P[%d] - P[%d] -> %.15lf\n", r2.index1, r2.index2, r2.Distanza);
printf("\nTempo impiegato per la ricerca -> %5.5f secondi\n", (double)(c_end - c_start) / CLOCKS_PER_SEC);

free(myArray.m);

return 0;
}

[gugoXX]

Cosa rappresenta quell'uno dopo il diecimila nella prima riga del file? :confused:

E' un centomila... Comunque il primo uno è il secondo parametro, il lato del cubo di contenimento, quindi è il valore massimo che puo' assumere la coordinata di uno qualsiasi dei punti.
Non sara' sempre 1 pertanto, potrebbe essere anche 3.14 volendo.
Ma l'avatar cos'è? Un nuovo tipo di faccia dopo l'ennesima tortura medioevale?

[Vincenzo1968]

E' un centomila... Comunque il primo uno è il secondo parametro, il lato del cubo di contenimento, quindi è il valore massimo che puo' assumere la coordinata di uno qualsiasi dei punti.
Non sara' sempre 1 pertanto, potrebbe essere anche 3.14 volendo.

Grazie mille :)


Ma l'avatar cos'è? Un nuovo tipo di faccia dopo l'ennesima tortura medioevale?
:D

No, è il bassista dei mitici Kiss:

http://it.youtube.com/watch?v=HnqUAbGMjLI

http://it.youtube.com/watch?v=DWLpbcgc814&feature=related

[VICIUS]

In sostanza per il primo punto si devono controllare tutte le coppie di punti presenti nel file e calcolare la distanza? Se è così è piuttosto facile da scrivere anche se se ho il sospetto che ruby impiegherà una vita per completare il file... :muro:

[gugoXX]

In sostanza per il primo punto si devono controllare tutte le coppie di punti presenti nel file e calcolare la distanza? Se è così è piuttosto facile da scrivere anche se se ho il sospetto che ruby impiegherà una vita per completare il file... :muro:

Certo, diciamo che questa è la soluzione di base di riferimento.
Si spera poi di trovare algoritmi migliori.

[VICIUS]

Certo, diciamo che questa è la soluzione di base di riferimento.
Si spera poi di trovare algoritmi migliori.
Ottimo. Ora sto provando a scrivere questa così ho qualcosa di corretto su cui basarmi.

[VICIUS]

Per ora sono arrivato a questo e come pensavo è di una lentezza disarmante.
Infinity = 1.0/0.0

class Point
attr_reader :x, :y, :z

def initialize(x, y, z)
@x, @y, @z = x.to_f, y.to_f, z.to_f
end

def distance_from(other)
# manca la radice
(@x - other.x) ** 2 +
(@y - other.y) ** 2 +
(@z - other.z) ** 2
end

def to_s
"(#{@x},#{@y},#{@z})"
end
end


class Result
attr_reader :a, :b, :dist

def initialize(a, b, dist)
@a, @b, @dist = a, b, dist
end

def to_s
"#{a} -> #{b} = #{Math.sqrt(dist)}"
end
end


def load_points_from_file(file_name)
list = Array.new
File.open(file_name, 'r') do |file|
file.gets # ignora la prima riga
while line = file.gets
list.push Point.new(*line.split)
end
end

list
end


def contest_9(list)
min, max = Result.new(0, 0, Infinity), Result.new(0, 0, -Infinity)
list.combination(2).each do |a,b|
dist = a.distance_from b
min = Result.new(a, b, dist) if dist < min.dist
max = Result.new(a, b, dist) if dist > max.dist
end

[min, max]
end

list = load_points_from_file('Coords.txt')
puts contest_9(list)

Con soli 2000 punti impiega circa 7 secondi quindi non chiedetemi di provare con quello completo che non ho alcuna intenzione di fondere il pc :)

[rеpne scasb]

■

[Vincenzo1968]

Accludo la mia versione:
...
Tempo per la creazione del software: meno di 30 minuti.
Algoritmo utilzzato: non e' la forza bruta.
Commento sul codice: e' scritto nel solito modo => incomprensibile :D
Commento sul contest: mi sarei aspettata qualcosa di meglio della forza bruta, cionci dove sei?

Io ho postato per primo quello a forza bruta in modo da seguire i criteri per la valutazione di cui si era parlato(e, cioè, prendere il rapporto fra il tempo della soluzione brute-force e quello della soluzione ottimizzata).


Addendum: a richesta, se interessa, accludero' una spiegazione dettagliata dell'algoritmo.

A me interessa. :)


P.S. Non mi pare ci siano bugs, ma nel caso ne troviate avvertitemi, grazie.

L'ho provato e funziona perfettamente:

http://www.guidealgoritmi.it/images/ImgForums/Contest08ASoluzRepne.jpg

[rеpne scasb]

■

[Vincenzo1968]

Per la mia soluzione ottimizzata ho ordinato codesto libro:

http://www.amazon.com/Computational-Geometry-Cambridge-Theoretical-Computer/dp/0521649765

:bimbo:

[rеpne scasb]

■

[Vincenzo1968]

E vabbe'...il tempo che ti arriva il libro, lo leggi, lo assimili, pensi ad una soluzione ottimizzata, la sviluppi....e' arrivato sicuramente Natale! :D

Qualcosa in tempi piu' rapidi?

eh lo so. :D
Dove posso documentarmi? Qualche link? C'entrano qualcosa i diagrammi di Voronoi? E i KD-Tree?

[rеpne scasb]

■

[Vincenzo1968]

Non saprei dirti, a naso mi pare di no. La mia attuale soluzione "ottimizzata", non e' nulla di particolare, e non servono libri o documentazione "particolare" per implementarla; non serve alcuna conoscenza matematica specifica.

La mia seconda soluzione ottimizzata (che non ho ancora pubblicato), sfrutta invece una proprieta' delle funzioni monotone crescenti. In questo caso sono necessarie alcune nozioni di matematica elementare.

C'entrano comunque gli alberi binari, giusto? (lo costruisci e successivamente lo linearizzi con Bin2Nor; è corretto?).

[rеpne scasb]

■

[Vincenzo1968]

Io provo coi kD-Trees utilizzando questa libreria:

http://code.google.com/p/kdtree/

:boxe:

[gugoXX]

Per ora sono qui.

ResMin : 70ms Dist: 0.000366913641159708 Indexes: 22503-72620
ResMax : 1546ms Dist: 1.70445415447397 Indexes: 72832-48483
1616ms

[Vincenzo1968]

Per ora sono qui.

ResMin : 70ms Dist: 0.000366913641159708 Indexes: 22503-72620
ResMax : 1546ms Dist: 1.70445415447397 Indexes: 72832-48483
1616ms


Io invece sono qui:

:lamer:

[rеpne scasb]

■

[gugoXX]

Non so, non sto utilizzando alberi di ricerca (Mi e' sembrato di aver capito che li stai utilizzando).
Forse riesco a tagliare ancora qualcosa, ma ci pensero' con calma.

Si', per la ricerca del massimo sto utilizzando un po' di parallelismo, ma non sono molto soddisfatto della soluzione per ora.

[rеpne scasb]

■

[Vincenzo1968]

Ohé,

sto provando con quelli che Knuth(The Art of Computer Programming (http://www-cs-faculty.stanford.edu/~knuth/taocp.html), Vol 3°, pag. 563) chiama post-office trees.

:lamer:

[Vincenzo1968]

Avrei bisogno di sapere come si fa a confrontare due punti.

Mi spiego:

data la seguente struttura:


typedef struct tagPunti
{
double A;
double B;
double C;
} Punti;


Dovrei scrivere una funzione


int Compare(Punti a, Punti b) // Funzione di confronto


che mi restituisca -1 se il punto a è minore del punto b; 1 se a > b e 0 se a == b.

Come si implementa una funzione di confronto tra punti con coordinate tridimensionali?

Grazie :)

[^TiGeRShArK^]

Avrei bisogno di sapere come si fa a confrontare due punti.

Mi spiego:

data la seguente struttura:


typedef struct tagPunti
{
double A;
double B;
double C;
} Punti;


Dovrei scrivere una funzione


int Compare(Punti a, Punti b) // Funzione di confronto


che mi restituisca -1 se il punto a è minore del punto b; 1 se a > b e 0 se a == b.

Come si implementa una funzione di confronto tra punti con coordinate tridimensionali?

Grazie :)

dovresti confrontare contemporaneamente le tre proiezioni sugli assi x, y, e z (che non sono altro che le coordinate del punto) credo....
EDIT:
e mi sa che devi anche stabilire una convenzione tua di ordinamento (per esempio assegnare una priorità diversa ai vari assi).
Solitamente si confrontano solo i moduli dei vettori definiti dai punti per ordinarli, per ordinare tutto il vettore devi stabilire tu se è + importante il modulo, o quale tra le due fasi....
O, come dicevo prima, un particolare asse...

[Vincenzo1968]

dovresti confrontare contemporaneamente le tre proiezioni sugli assi x, y, e z (che non sono altro che le coordinate del punto) credo....

Grazie Tiger ;)

per il momento sono riuscito ad ottenere questo:

http://www.guidealgoritmi.it/images/ImgForums/Contest08Soluz01OptTemp.jpg

I tempi sono buoni ma la funzione di confronto(evidentemente sbagliata) mi fa ottenere risultati errati :cry:

Non è che puoi postarmi uno pseudo codice?

:bimbo:

[Tommo]

Uff... vorrei provare anche io, ma ho strani problemi di null pointers :stordita:

Cmq volevo chiedere, a che serve usare un'albero di ricerca quando i dati sono ben conosciuti e puntiformi?

Non basta suddividere lo spazio in cellette messe in un array normale?



Cell* myCell = cells[ pos.x + side.x*pos.y + side.x*side.y*pos.z ];

Particle** nearParticles = myCell->getContainedParticles();
size_t nearParticlesNb = myCell->getContainedParticlesNb();

//confronta la distanza minima solo con le celle più vicine


Il vantaggio è che al crescere del numero di celle, trovare quella che contiene la particella costa esattamente uguale.
In più, mediamente le particelle in una cella sono un numero costante, quindi se celle/particelle si mantiene costante, il tutto pesa sempre uguale. Quasi.
Però: non funge nel caso che le due particelle più vicine siano a cavallo di due celle. :stordita:

Ha senso?

[^TiGeRShArK^]

Grazie Tiger ;)

per il momento sono riuscito ad ottenere queso:

http://www.guidealgoritmi.it/images/ImgForums/Contest08Soluz01OptTemp.jpg

I tempi sono buoni ma la funzione di confronto(evidentemente sbagliata) mi fa ottenere risultati errati :cry:

Non è che puoi postarmi uno pseudo codice?

:bimbo:

eh beh..
non è che è errata la funzione di confronto è che non penso esista una funzione di confronto univocamente definita per dei vettori in uno spazio euclideo.. :stordita:
Tanto per farti un esempio se confronti semplicemente le tre coordinate senza dare una priorità come ti comporti nel caso in cui i due punti siano (1,5,2) e (2,3,1)?
la coordinata x è minore nel primo punto, ma la y e la z sono maggiori...
dunque se decidi che a parità di x vai a confrontare la y e quindi a parità di entrambi confronti la z puoi fare un ordinamento.
in pratica in pseudo-codice sarebbe qualcosa del genere:

def compare(p1, p2) {
if (p1.x < p2.x || p1.y < p2.y || p1.z < p2.z)
return -1
else if (p1.x == p2.x && p1.y == p2.y && p1.z == p2.z)
return 0
else
return 1

in questo modo tutti i punti in cui almeno una delle coordinate è minore della corrispettiva coordinata del punto da confrontare sono minori di esso, e sono maggiori se tutte e tre le coordinate sono maggiori.
Ma si deve vedere se a te va bene questo tipo di ordinamento o se ne devi utilizzare qualche altro.. :stordita:

[Vincenzo1968]

Posto il codice:


#include <stdio.h>
#include <math.h>
#include <time.h>
#include <malloc.h>

#define BUFFER_SIZE 4096

int g_pass = 2;

//double distanzaMin = -1;
//double distanzaMax = -1;

typedef struct tagPunti
{
int index;
double A;
double B;
double C;
} Punti;

typedef struct tagRes
{
Punti P1;
Punti P2;
double Distanza;
int index1;
int index2;
} Res;

typedef struct tagArray
{
int dimensione;
double dmax;
Punti *m;
} Array;

typedef enum tagStati
{
S_ERROR = -1, S0 = 0, S1, S2, S3, S4, S5, S6, S7, S8, S9
} Stati;


// La macro seguente restituisce il numero degli elementi di un array.
// Per esempio, data la seguente dichiarazione:
//
// int array[] = {8, 5, 13, 2, 1};
//
// ARRAY_SIZE(array) restituisce 5

#define ARRAY_SIZE(Array) (sizeof(Array) / sizeof((Array)[0]))

typedef Punti Item;
//typedef double Item;

int (*pfnCompare)(Item, Item) = NULL; // Puntatore alla funzione di confronto

int Compare(Item a, Item b) // Funzione di confronto
{
// La funzione di confronto deve restituire un valore
// minore di zero se a < b,
// maggiore di zero se a > b
// e zero se a e b sono uguali

double x, y;

if ( g_pass == 0 )
{
x = a.A;
y = b.A;
g_pass = 1;
}
else if ( g_pass == 1 )
{
x = a.B;
y = b.B;
g_pass = 2;
}
else if ( g_pass == 2 )
{
x = a.C;
y = b.C;
g_pass = 0;
}

if ( x < y )
return -1;
else if ( x > y )
return 1;
else
return 0;
}

Stati DFA(const char *szFileName, Array *pArray)
{
Stati stato = S0;
FILE *fp;
unsigned char buffer[BUFFER_SIZE];
int numblocks;
int numread;
unsigned char c;
int k, x, j;
int riga, colonna;
char szNum[256];
double num;

unsigned char byteCR = 0xD; // Carriage Return
unsigned char byteLF = 0xA; // Line Feed

fp = fopen(szFileName, "rb");
if ( fp == NULL )
{
printf("Errore nell'apertura del file %s\n", szFileName);
return S_ERROR;
}

if ( fseek(fp, 0, SEEK_END) )
return 0;

numblocks = ftell(fp)/BUFFER_SIZE;
if ( numblocks == 0 )
{
numblocks = 1;
}
else
{
if ( ftell(fp) % BUFFER_SIZE != 0 )
numblocks++;
}

fseek(fp, 0, SEEK_SET);
numread = fread(buffer, 1, BUFFER_SIZE, fp);
if (numread == 0)
{
fclose(fp);
printf("Errore 1 nella lettura del file %s\n", szFileName);
return S_ERROR;
}

pArray->dimensione = 0;
k = 0;
c = *(buffer + k++);
while ( 1 )
{
if ( c >= '0' && c <= '9' )
{
pArray->dimensione = pArray->dimensione * 10 + c - '0';
}
else if ( c == ' ' || c == '\t' )
{
j = 0;
while ( c != byteLF )
{
c = *(buffer + k++);
if ( c >= '0' && c <= '9' )
szNum[j++] = c;
}
szNum[j] = '\0';
pArray->dmax = atof(szNum);
break;
}
else if ( c == '\n' )
{
break;
}

c = *(buffer + k++);
}

pArray->m = (Punti*)malloc(pArray->dimensione*sizeof(Punti));
if ( !pArray->m )
{
printf("Memoria insufficiente.\n");
fclose(fp);
return S_ERROR;
}

riga = colonna = 0;
num = 0;
x = j = 0;
while ( x < numblocks )
{
c = *(buffer + k++);

if ( c == byteLF )
{
pArray->m[riga].index = riga;
riga++;
colonna = 0;
}

switch (stato)
{
case S0:
j = 0;
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S2;
}
else if ( c == '.' )
{
szNum[j++] = c;
stato = S3;
}
else if ( c == '-' || c == '+' )
{
szNum[j++] = c;
stato = S1;
}
else if ( c == ' ' || c == '\t' )
{
while ( c == ' ' || c == '\t' )
{
c = *(buffer + k++);
if ( k >= numread )
break;
}
k--;
}
else if ( c == byteCR || c == byteLF )
{
// nada
}
else
{
printf("Errore: stato S0; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S1:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S2;
}
else if ( c == '.' )
{
szNum[j++] = c;
stato = S3;
}
else
{
printf("Errore: stato S1; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S2:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
}
else if ( c == '.' )
{
szNum[j++] = c;
stato = S4;
}
else if ( c == 'E' || c == 'e' )
{
szNum[j++] = c;
stato = S5;
}
else if ( c == ' ' || c == '\t' || c == byteCR || c == byteLF )
{
szNum[j] = '\0';
if ( colonna == 0 )
pArray->m[riga].A = atof(szNum);
else if ( colonna == 1 )
pArray->m[riga].B = atof(szNum);
else if ( colonna == 2 )
pArray->m[riga].C = atof(szNum);
colonna++;
stato = S0;
}
else
{
printf("Errore: stato S2; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S3:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S6;
}
else
{
printf("Errore: stato S3; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S4:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S6;
}
else if ( c == 'E' || c == 'e' )
{
szNum[j++] = c;
stato = S7;
}
else if ( c == ' ' || c == '\t' || c == byteCR || c == byteLF )
{
szNum[j] = '\0';
if ( colonna == 0 )
pArray->m[riga].A = atof(szNum);
else if ( colonna == 1 )
pArray->m[riga].B = atof(szNum);
else if ( colonna == 2 )
pArray->m[riga].C = atof(szNum);
colonna++;
stato = S0;
}
else
{
printf("Errore: stato S4; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S5:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S9;
}
else if ( c == '-' || c == '+' )
{
szNum[j++] = c;
stato = S8;
}
else
{
printf("Errore: stato S5; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S6:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S6;
}
else if ( c == 'E' || c == 'e' )
{
szNum[j++] = c;
stato = S7;
}
else if ( c == ' ' || c == '\t' || c == byteCR || c == byteLF )
{
szNum[j] = '\0';
if ( colonna == 0 )
pArray->m[riga].A = atof(szNum);
else if ( colonna == 1 )
pArray->m[riga].B = atof(szNum);
else if ( colonna == 2 )
pArray->m[riga].C = atof(szNum);
colonna++;
stato = S0;
}
else
{
printf("Errore: stato S6; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S7:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S9;
}
else if ( c == '-' || c == '+' )
{
szNum[j++] = c;
stato = S8;
}
else
{
printf("Errore: stato S7; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S8:
case S9:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S9;
}
else if ( c == ' ' || c == '\t' || c == byteCR || c == byteLF )
{
szNum[j] = '\0';
if ( colonna == 0 )
pArray->m[riga].A = atof(szNum);
else if ( colonna == 1 )
pArray->m[riga].B = atof(szNum);
else if ( colonna == 2 )
pArray->m[riga].C = atof(szNum);
colonna++;
stato = S0;
}
else
{
printf("Errore: stato S9; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;
}

if ( k >= numread )
{
numread = fread(buffer, 1, BUFFER_SIZE, fp);
k = 0;
x++;
}
}

fclose(fp);

return stato;
}

void Calcola(Punti P[], int dim, Res *pMin, Res *pMax)
{
int x;
double distanza, distanzaMin, distanzaMax;
double diffA, diffB, diffC;

distanzaMin = distanzaMax = -1;

for ( x = 1; x < dim; x++ )
{
diffA = P[x-1].A - P[x].A;
diffB = P[x-1].B - P[x].B;
diffC = P[x-1].C - P[x].C;
distanza = diffA*diffA + diffB*diffB + diffC*diffC;
if ( distanza < distanzaMin || distanzaMin == -1 )
{
distanzaMin = distanza;
pMin->P1 = P[x-1];
pMin->P2 = P[x];
pMin->index1 = P[x-1].index;
pMin->index2 = P[x].index;
}
if ( distanza > distanzaMax )
{
distanzaMax = distanza;
pMax->P1 = P[x-1];
pMax->P2 = P[x];
pMax->index1 = P[x-1].index;
pMax->index2 = P[x].index;
}
}

pMin->Distanza = sqrt( distanzaMin );
pMax->Distanza = sqrt( distanzaMax );
}

void QuickSortNonRecursive(Item *a, int start, int stop)
{
int i, s = 0, stack[64];

int up = start, down = stop - 1;
Item part = a[stop];
Item tmp;

stack[s++] = start;
stack[s++] = stop;

while (s > 0)
{
stop = stack[--s];
start = stack[--s];

if (start >= stop)
continue;

//INIZIO PARTIZIONAMENTO

up = start;
down = stop - 1;
part = a[stop];

if (stop <= start)
{
i = start;
continue;
}

while (1)
{
while ( (*pfnCompare)(a[up], part) < 0 )
up++;

while ( (*pfnCompare)(part, a[down]) < 0 && (up < down) )
down--;

if (up >= down)
break;

tmp = a[down];
a[down] = a[up];
a[up] = tmp;

up++;
down--;
}

tmp = a[stop];
a[stop] = a[up];
a[up] = tmp;

i = up;

// FINE PARTIZIONAMENTO

if (i - start > stop - i)
{
stack[s++] = start;
stack[s++] = i - 1;
stack[s++] = i + 1;
stack[s++] = stop;
}
else
{
stack[s++] = i + 1;
stack[s++] = stop;
stack[s++] = start;
stack[s++] = i - 1;
}
}
}

int main()
{
Stati stato;
Array myArray;
Res r1, r2;

clock_t c_start, c_end;

char *szFileName = "C:\\Scaricamenti\\Temp\\Gugo\\Contest 08 - Compattamento Barionico\\Coords\\Coords.dat";

pfnCompare = Compare;

c_start = clock();
stato = DFA(szFileName, &myArray);
c_end = clock();
printf("\nTempo impiegato per caricamento dati -> %5.5f secondi\n", (double)(c_end - c_start) / CLOCKS_PER_SEC);

c_start = clock();
QuickSortNonRecursive(myArray.m, 0, myArray.dimensione - 1);
Calcola(myArray.m, myArray.dimensione, &r1, &r2);
c_end = clock();

printf("\nDistanza Min: P[%d] - P[%d] -> %.15lf\n", r1.index1, r1.index2, r1.Distanza);
printf("\nDistanza Max: P[%d] - P[%d] -> %.15lf\n", r2.index1, r2.index2, r2.Distanza);
printf("\nTempo impiegato per la ricerca -> %5.5f secondi\n", (double)(c_end - c_start) / CLOCKS_PER_SEC);

free(myArray.m);

return 0;
}


Per ottenere risultati corretti bisogna modificare la funzione Compare

:bimbo:

[Vincenzo1968]

eh beh..
non è che è errata la funzione di confronto è che non penso esista una funzione di confronto univocamente definita per dei vettori in uno spazio euclideo.. :stordita:
Tanto per farti un esempio se confronti semplicemente le tre coordinate senza dare una priorità come ti comporti nel caso in cui i due punti siano (1,5,2) e (2,3,1)?
la coordinata x è minore nel primo punto, ma la y e la z sono maggiori...
dunque se decidi che a parità di x vai a confrontare la y e quindi a parità di entrambi confronti la z puoi fare un ordinamento.
in pratica in pseudo-codice sarebbe qualcosa del genere:

def compare(p1, p2) {
if (p1.x < p2.x || p1.y < p2.y || p1.z < p2.z)
return -1
else if (p1.x == p2.x && p1.y == p2.y && p1.z == p2.z)
return 0
else
return 1

in questo modo tutti i punti in cui almeno una delle coordinate è minore della corrispettiva coordinata del punto da confrontare sono minori di esso, e sono maggiori se tutte e tre le coordinate sono maggiori.
Ma si deve vedere se a te va bene questo tipo di ordinamento o se ne devi utilizzare qualche altro.. :stordita:

Io basavo il mio ragionamento sul caso unidimensionale.
Mi spiego: se ho un array non ordinato di double(che rappresenta un insieme di punti su una retta), per ricercare i due punti a distanza minima, debbo confrontare tutte le possibili coppie di punti(complessità = O(n^2).
Se invece ordino l'array, per trovare la coppia di punti a distanza minima mi basta scorrere l'array confrontando gli elementi adiacenti(complessità = O(n)).

Ho detto minchiate?

Se il ragionamento sul caso unidimensionale è corretto, è possibile applicarlo anche al caso tridimensionale?

[gugoXX]

Io basavo il mio ragionamento sul caso unidimensionale.
Mi spiego: se ho un array non ordinato di double(che rappresenta un insieme di punti su una retta), per ricercare i due punti a distanza minima, debbo confrontare tutte le possibili coppie di punti(complessità = O(n^2).
Se invece ordino l'array, per trovare la coppia di punti a distanza minima mi basta scorrere l'array confrontando gli elementi adiacenti(complessità = O(n)).

Ho detto minchiate?

Se il ragionamento sul caso unidimensionale è corretto, è possibile applicarlo anche al caso tridimensionale?

Purtroppo no, non e' applicabile neppure al caso bidimensionale.
Se la distanza PQ e' 5, e la distanza PR e' 6, non e' detto che QR sia solo 1.
Non sono neppure mantenute le distanze reciproche. Magari Q ed R distano tra loro piu' di quanto distavano da P (es: Triangolo ottusangolo)

Quindi non puoi ordinare i punti di un piano lungo una retta, pensando di calcolare le distanze reciproche solo ragionando sulla posizione del punto lungo la retta.
Occorre un'altra strategia.

[VICIUS]

Purtroppo no, non e' applicabile neppure al caso bidimensionale.
Se la distanza PQ e' 5, e la distanza PR e' 6, non e' detto che QR sia solo 1.
Non sono neppure mantenute le distanze reciproche. Magari Q ed R distano tra loro piu' di quanto distavano da P (es: Triangolo ottusangolo)

Quindi non puoi ordinare i punti di un piano lungo una retta, pensando di calcolare le distanze reciproche solo ragionando sulla posizione del punto lungo la retta.
Occorre un'altra strategia.
Immagino sia dovuto alla perdita dell'ordinamento nel passaggio da R a C giusto?

[Vincenzo1968]

Purtroppo no, non e' applicabile neppure al caso bidimensionale.
Se la distanza PQ e' 5, e la distanza PR e' 6, non e' detto che QR sia solo 1.
Non sono neppure mantenute le distanze reciproche. Magari Q ed R distano tra loro piu' di quanto distavano da P (es: Triangolo ottusangolo)

Quindi non puoi ordinare i punti di un piano lungo una retta, pensando di calcolare le distanze reciproche solo ragionando sulla posizione del punto lungo la retta.
Occorre un'altra strategia.

Ok, grazie mille ;)

Mò provo con i kD-Trees e buona notte ai suonatori.

:bimbo:

[^TiGeRShArK^]

Immagino sia dovuto alla perdita dell'ordinamento nel passaggio da R a C giusto?

più che altro mi sa che R è isomorfo solo con una delle tre basi di R^3 che è il nostro spazio euclideo...
Quindi non è possibile stabilire un isomorfismo tra l'intero spazio vettoriale e il solo R e mi sa che così si perde l'ordinabilità assoluta.....
ma prendete con le pinze quello che ho detto che l'algebra lineare mi è sempre stata sulle bolls.. :stordita:

[gugoXX]

Immagino sia dovuto alla perdita dell'ordinamento nel passaggio da R a C giusto?

Diciamo che e' simile al passaggio da C ad R.
Non puoi trovare una relazione d'ordine tra i numeri complessi, cosi' come non la puoi trovare tra i punti del piano.

Non puoi dire che il punto P sia maggiore del punto Q o viceversa.

PS: Ho pensato un algoritmo greedy* che dovrebbe essere decisamente veloce. Stasera lo posto

greedy*: Algoritmo tipicamente molto veloce (o semplice) che riesce a trovare una soluzione soddisfacente ma che non e' detto che sia la soluzione ottima.
http://it.wikipedia.org/wiki/Algoritmo_greedy

Un algoritmo greedy è un algoritmo che cerca di ottenere una soluzione ottima da un punto di vista globale attraverso la scelta della soluzione più golosa (aggressiva o avida, a seconda della traduzione preferita del termine greedy dall'inglese) ad ogni passo locale. Questa tecnica consente, dove applicabile (infatti non sempre si arriva ad una soluzione ottima), di trovare soluzioni ottimali per determinati problemi in un tempo polinomiale (cfr. Problemi NP-Completi, cioè problemi di soluzione non deterministica polinomiale).

Vedremo il risultato.

[marco.r]

Uhm, contest poco utile secondo me perche' c'e' sia un sacco di teoria che di codice gia' fatto per trovare i punti piu' vicini.
Comunque, visto che il bello dei contest e' che servono a poco e ci si diverte un sacco, posto anche io una mia soluzione custom, per il momento per la coppia piu' lontana.


#include <iostream>
#include <vector>
#include <fstream>
#include <set>
#include <iterator>
#include <cmath>
#include <algorithm>

using std::cout;
using std::endl;
using std::make_pair;
using std::set;
using std::pair;

struct point
{
point(){}
point(double _x, double _y, double _z ):x(_x),y(_y),z(_z){}
double x,y,z;
};

std::vector<point> data;
std::set<int> feasible_set;
point corners[8];

double max_distance = 0.0;
std::pair<int,int> max_pair;
double size;
int dim;



std::istream& operator >> (std::istream& in, point& p )
{
in >> p.x >> p.y >> p.z;
return in;
}

std::ostream& operator << (std::ostream& out, point& p )
{
out << p.x << ' ' << p.y << ' ' << p.z;
return out;
}


double distance(const point& p1, const point& p2 )
{
double dx,dy,dz;
dx = p1.x - p2.x;
dy = p1.y - p2.y;
dz = p1.z - p2.z;
return dx*dx + dy*dy + dz*dz;
}

pair<int,double> most_distant(const point& p)
{
double dist = 0;
int idx = -1;
for ( std::set<int>::const_iterator i = feasible_set.begin();
i != feasible_set.end();
++i )
{
double d = distance(p,data[*i]);
if ( d > dist )
{
dist = d;
idx = *i;
}
}
return make_pair(idx,dist);
}

bool feasible( const point& p )
{
for ( int i=0 ; i<8 ; ++i )
{
if ( distance(p,corners[i]) > max_distance )
return true;
}
return false;

}


void init( double size )
{
corners[0] = point( 0, 0, 0 );
corners[1] = point( 0, 0, 1 );
corners[2] = point( 0, 1, 0 );
corners[3] = point( 0, 1, 1 );
corners[4] = point( 1, 0, 0 );
corners[5] = point( 1, 0, 1 );
corners[6] = point( 1, 1, 0 );
corners[7] = point( 1, 1, 1 );

}

bool unfeasible( const int& x )
{
return !feasible(data[x]);
}


void remove_unfeasible()
{
for ( set<int>::const_iterator i = feasible_set.begin();
i != feasible_set.end();
++i )
{
if ( unfeasible(*i) )
{
feasible_set.erase( i );
}
}
}

void update_feasible_set( const point& p, int index )
{
if ( feasible(p) )
{
pair<int,double> x = most_distant(p);
const int& idx = x.first;
const double& d = x.second;
if ( d > max_distance )
{
max_distance = d;
max_pair = std::make_pair( index, idx );
}

remove_unfeasible();
feasible_set.insert(index);
}
}

void solve_max()
{
std::cout << "Looking for greatest distance" << std::endl;
feasible_set.insert( 0 );
feasible_set.insert( 1 );
max_distance = distance( data[0], data[1] );
max_pair = make_pair(0,1);
for ( int i=2 ; i<data.size() ; ++i )
{
update_feasible_set( data[i], i );
}
cout << "Max distance is " << sqrt(max_distance) << " [" << max_pair.first << ',' << max_pair.second << "]" << endl;
cout << data[max_pair.first] << endl;
cout << data[max_pair.second] << endl;
}

int main(int argc, char** argv)
{
const char* filename = argc > 1 ? argv[1] : "Coords.dat";

std::ifstream input(filename);
std::cout << "Reading data set ... " << std::endl;
input >> dim >> size;
init(size);

data.resize(dim);
std::copy( std::istream_iterator<point>(input),
std::istream_iterator<point>(),
data.begin() );

solve_max();

}

L'idea di base e' che i punti piu' distanti staranno al'langolo del cubo. Per cui an mano che scorro i punti letti elimino quelli che non possono contribuire alla soluzione finale perche' troppo in centro (ovvero distano dall'angolo del cubo piu' lontano meno della soluzione corrente migliore).
Non pensavo, ma anche con diverse centinaia di migliaia di valori, riesco a tenere il working set sotto la cinquantina di elementi, per cui anche un confronto "brute force" con questi risulta sufficiente.
Ho una idea custom anche per l'altro problema, stasera che avro' tempo spero di svilupparla un attimo.

[marco.r]

ops,
ovviamente c'e' un errore nella init :D, irrilevante per il dataset attuale cmq (che sta nel cubo unitario)

[Vincenzo1968]

...
PS: Ho pensato un algoritmo greedy* che dovrebbe essere decisamente veloce. Stasera lo posto

greedy*: Algoritmo tipicamente molto veloce (o semplice) che riesce a trovare una soluzione soddisfacente ma che non e' detto che sia la soluzione ottima.
http://it.wikipedia.org/wiki/Algoritmo_greedy

Vedremo il risultato.

Io ce l'ho un algoritmo greedy ma non è molto veloce(ci vogliono due giorni):

Tratto da Gli arancini di Montalbano di Andrea Camilleri:

Adelina ci metteva due jornate sane sane a pripararli. Ne sapeva, a memoria, la ricetta. Il giorno avanti si fa un aggrassato di vitellone e di maiale in parti uguali che deve còciri a foco lentissimo per ore e ore con cipolla, pummadoro, sedano, prezzemolo e basilico. Il giorno appresso si pripara un risotto, quello che chiamano alla milanisa(senza zaffirano, pi carità!), lo si versa sopra a una tavola, ci si impastano le ova e lo si fa rifriddare. Intanto si còcino i pisellini, si fa una besciamella, si riducono a pezzettini 'na poco di fette di salame e si fa tutta una composta con la carne aggrassata, triturata a mano con la mezzaluna(nenti frullatore, pi carità di Dio!). Il suco della carne s'ammisca col risotto. A questo punto si piglia tanticchia di risotto, s'assistema nel palmo d'una mano fatta a conca, ci si mette dentro quanto un cucchiaio di composta e si copre con dell'altro riso a formare una bella palla. Ogni palla la si fa rotolare nella farina, poi si passa nel bianco d'ovo e nel pane grattato. Doppo, tutti gli arancini s'infilano in una paredda d'oglio bollente e si fanno friggere fino a quando pigliano un colore d'oro vecchio. Si lasciano scolare sulla carta. E alla fine, ringraziannu u Signuruzzu, si mangiano!

[^TiGeRShArK^]

Uhm, contest poco utile secondo me perche' c'e' sia un sacco di teoria che di codice gia' fatto per trovare i punti piu' vicini.
Comunque, visto che il bello dei contest e' che servono a poco e ci si diverte un sacco, posto anche io una mia soluzione custom, per il momento per la coppia piu' lontana.


#include <iostream>
#include <vector>
#include <fstream>
#include <set>
#include <iterator>
#include <cmath>
#include <algorithm>

using std::cout;
using std::endl;
using std::make_pair;
using std::set;
using std::pair;

struct point
{
point(){}
point(double _x, double _y, double _z ):x(_x),y(_y),z(_z){}
double x,y,z;
};

std::vector<point> data;
std::set<int> feasible_set;
point corners[8];

double max_distance = 0.0;
std::pair<int,int> max_pair;
double size;
int dim;



std::istream& operator >> (std::istream& in, point& p )
{
in >> p.x >> p.y >> p.z;
return in;
}

std::ostream& operator << (std::ostream& out, point& p )
{
out << p.x << ' ' << p.y << ' ' << p.z;
return out;
}


double distance(const point& p1, const point& p2 )
{
double dx,dy,dz;
dx = p1.x - p2.x;
dy = p1.y - p2.y;
dz = p1.z - p2.z;
return dx*dx + dy*dy + dz*dz;
}

pair<int,double> most_distant(const point& p)
{
double dist = 0;
int idx = -1;
for ( std::set<int>::const_iterator i = feasible_set.begin();
i != feasible_set.end();
++i )
{
double d = distance(p,data[*i]);
if ( d > dist )
{
dist = d;
idx = *i;
}
}
return make_pair(idx,dist);
}

bool feasible( const point& p )
{
for ( int i=0 ; i<8 ; ++i )
{
if ( distance(p,corners[i]) > max_distance )
return true;
}
return false;

}


void init( double size )
{
corners[0] = point( 0, 0, 0 );
corners[1] = point( 0, 0, 1 );
corners[2] = point( 0, 1, 0 );
corners[3] = point( 0, 1, 1 );
corners[4] = point( 1, 0, 0 );
corners[5] = point( 1, 0, 1 );
corners[6] = point( 1, 1, 0 );
corners[7] = point( 1, 1, 1 );

}

bool unfeasible( const int& x )
{
return !feasible(data[x]);
}


void remove_unfeasible()
{
for ( set<int>::const_iterator i = feasible_set.begin();
i != feasible_set.end();
++i )
{
if ( unfeasible(*i) )
{
feasible_set.erase( i );
}
}
}

void update_feasible_set( const point& p, int index )
{
if ( feasible(p) )
{
pair<int,double> x = most_distant(p);
const int& idx = x.first;
const double& d = x.second;
if ( d > max_distance )
{
max_distance = d;
max_pair = std::make_pair( index, idx );
}

remove_unfeasible();
feasible_set.insert(index);
}
}

void solve_max()
{
std::cout << "Looking for greatest distance" << std::endl;
feasible_set.insert( 0 );
feasible_set.insert( 1 );
max_distance = distance( data[0], data[1] );
max_pair = make_pair(0,1);
for ( int i=2 ; i<data.size() ; ++i )
{
update_feasible_set( data[i], i );
}
cout << "Max distance is " << sqrt(max_distance) << " [" << max_pair.first << ',' << max_pair.second << "]" << endl;
cout << data[max_pair.first] << endl;
cout << data[max_pair.second] << endl;
}

int main(int argc, char** argv)
{
const char* filename = argc > 1 ? argv[1] : "Coords.dat";

std::ifstream input(filename);
std::cout << "Reading data set ... " << std::endl;
input >> dim >> size;
init(size);

data.resize(dim);
std::copy( std::istream_iterator<point>(input),
std::istream_iterator<point>(),
data.begin() );

solve_max();

}

L'idea di base e' che i punti piu' distanti staranno al'langolo del cubo. Per cui an mano che scorro i punti letti elimino quelli che non possono contribuire alla soluzione finale perche' troppo in centro (ovvero distano dall'angolo del cubo piu' lontano meno della soluzione corrente migliore).
Non pensavo, ma anche con diverse centinaia di migliaia di valori, riesco a tenere il working set sotto la cinquantina di elementi, per cui anche un confronto "brute force" con questi risulta sufficiente.
Ho una idea custom anche per l'altro problema, stasera che avro' tempo spero di svilupparla un attimo.
ecco, bravo...
hai scritto quello su cui stavo lavorando io :fagiano:
io però lo stavo facendo con le coordinate sferiche :stordita:

[VICIUS]

più che altro mi sa che R è isomorfo solo con una delle tre basi di R^3 che è il nostro spazio euclideo [...]

Diciamo che e' simile al passaggio da C ad R.
Non puoi trovare una relazione d'ordine tra i numeri complessi, cosi' come non la puoi trovare tra i punti del piano. [...]
Capito. Vado a rileggermi gli appunti di algebra 1... se li trovo ancora :asd:

[rеpne scasb]

■

[marco.r]

Ecco la mia versione.
Fondamentalmente sfrutto il fatto che i valori sono distribuiti uniformemente e che ho uno scatolo di lato prefissato per suddividere lo spazio e limitare i confronti "brute force" tra sottovolumi di piccola dimensione.
Scegliendo la dimensione in modo che il numero medio di punti sia costante al crescere del numero di punti totale si riesce a tenere sotto controllo la complessita'.
Soluzione un po' becera ma che comunque funziona. I tempi sono all'incirca quelli della soluzione di repne.

#include <iostream>
#include <vector>
#include <fstream>
#include <set>
#include <iterator>
#include <cmath>
#include <algorithm>
#include <cassert>

using std::cout;
using std::endl;
using std::make_pair;
using std::set;
using std::pair;
using std::vector;

clock_t start;

void start_clock()
{
start = clock();
}

double stop_clock()
{
clock_t new_start = clock();
float result = ((double)(new_start - start))/CLOCKS_PER_SEC;
start = new_start;
return result;
}

struct point
{
point(){}
point(double _x, double _y, double _z ):x(_x),y(_y),z(_z){}
double x,y,z;
};

std::vector<point> data;
std::set<int> feasible_set;
point corners[8];

double max_distance = 0.0;
std::pair<int,int> max_pair;
double size;
int dim;



std::istream& operator >> (std::istream& in, point& p )
{
in >> p.x >> p.y >> p.z;
return in;
}

std::ostream& operator << (std::ostream& out, const point& p )
{
out << p.x << ' ' << p.y << ' ' << p.z;
return out;
}


double distance(const point& p1, const point& p2 )
{
double dx,dy,dz;
dx = p1.x - p2.x;
dy = p1.y - p2.y;
dz = p1.z - p2.z;
return dx*dx + dy*dy + dz*dz;
}

pair<int,double> most_distant(const point& p)
{
double dist = 0;
int idx = -1;
for ( std::set<int>::const_iterator i = feasible_set.begin();
i != feasible_set.end();
++i )
{
double d = distance(p,data[*i]);
if ( d > dist )
{
dist = d;
idx = *i;
}
}
return make_pair(idx,dist);
}

bool feasible( const point& p )
{
for ( int i=0 ; i<8 ; ++i )
{
if ( distance(p,corners[i]) > max_distance )
return true;
}
return false;

}


void init( double size )
{
corners[0] = point( 0, 0, 0 );
corners[1] = point( 0, 0, 1 );
corners[2] = point( 0, 1, 0 );
corners[3] = point( 0, 1, 1 );
corners[4] = point( 1, 0, 0 );
corners[5] = point( 1, 0, 1 );
corners[6] = point( 1, 1, 0 );
corners[7] = point( 1, 1, 1 );

}

bool unfeasible( const int& x )
{
return !feasible(data[x]);
}


void remove_unfeasible()
{
for ( set<int>::const_iterator i = feasible_set.begin();
i != feasible_set.end();
++i )
{
if ( unfeasible(*i) )
{
feasible_set.erase( i );
}
}
}

void update_feasible_set( const point& p, int index )
{
if ( feasible(p) )
{
pair<int,double> x = most_distant(p);
const int& idx = x.first;
const double& d = x.second;
if ( d > max_distance )
{
max_distance = d;
max_pair = std::make_pair( index, idx );
}

remove_unfeasible();
feasible_set.insert(index);
}
}

void solve_max()
{
feasible_set.insert( 0 );
feasible_set.insert( 1 );
max_distance = distance( data[0], data[1] );
max_pair = make_pair(0,1);
for ( int i=2 ; i<data.size() ; ++i )
{
update_feasible_set( data[i], i );
}
cout << "Maximum distance is " << sqrt(max_distance) << " between [" << max_pair.first << ',' << max_pair.second << "]" << endl;
}


int nblocks( int points )
{
// voglio avere mediamente 3 punti per blocco

return int( pow( points/3.0, 1.0/3.0 ) )+1;
}

float blocksize( int points, float size )
{
return size/nblocks(points);
}

int block_index(int blocknum , int xidx,int yidx, int zidx)
{
const int& result = xidx*blocknum*blocknum + yidx*blocknum + zidx;
return result;
}
int point_block( float size, float blocksize, int blocknum, const point& p )
{
int xidx,yidx,zidx;
xidx = int( p.x * size/blocksize );
yidx = int( p.y * size/blocksize );
zidx = int( p.z * size/blocksize );
return block_index(blocknum,xidx,yidx,zidx);
}


void throw_in_buckets( vector< set<int> >& buckets, const vector<point>& points,
float size, float blocksize, int blocknum )
{
for ( int i=0 ; i<points.size() ; ++i )
{
int n = point_block( size, blocksize, blocknum, points[i] );
buckets[n].insert(i);
}
}

void find_closest_pair( const vector< set<int> >& buckets, int index1, int index2,
float& min_distance, int& first_index, int& second_index )
{
const set<int>& from = buckets[index1];
const set<int>& to = buckets[index2];

for ( set<int>::const_iterator i = from.begin();
i != from.end();
++i )
{
for ( set<int>::const_iterator j = to.begin();
j != to.end();
++j )
{
if ( *i == *j )
continue;
float d = distance( data[*i], data[*j] );
if ( d < min_distance )
{
first_index = *i;
second_index = *j;
min_distance = d;
}
}
}
}



void solve_min()
{
// totale confronti : N^2
// se divido in K blocchi con in mesia (N/K) valori:
// K * (N/K)^2 confronti = N^2/K
// prendo K = N/100 = N

int blocknum = nblocks(dim);
float bsize = blocksize(dim,size);
vector< set<int> > buckets;
buckets.resize( blocknum*blocknum*blocknum );

throw_in_buckets( buckets, data, size, bsize, blocknum );

const int deltas[8][3] = { {0,0,0},{0,0,1},{0,1,0},{0,1,1},
{1,0,0},{1,0,1},{1,1,0},{1,1,1}};
int first_index = 0;
int second_index = 1;
float min_distance = distance( data[0], data[1] );
for ( int i=0 ; i<blocknum-1 ; ++i )
{
for ( int j=0 ; j<blocknum-1; ++j )
{
for ( int k=0; k<blocknum-1; ++k )
{
for ( int n=0 ; n<8; ++n )
{
int index1 = block_index( blocknum,i,j,k );
int index2 = block_index( blocknum,
i+deltas[n][0],
j+deltas[n][1],
k+deltas[n][2] );
find_closest_pair( buckets, index1, index2, min_distance, first_index, second_index );
}
}
}
}
cout << "Minimum distance is " << sqrt(min_distance) << " between [" << first_index << ',' << second_index << "]" << endl;
}

int main(int argc, char** argv)
{
const char* filename = argc > 1 ? argv[1] : "Coords.dat";

std::ifstream input(filename);
float load_time,min_time,max_time;

start_clock();
input >> dim >> size;
init(size);

data.resize(dim);
std::copy( std::istream_iterator<point>(input),
std::istream_iterator<point>(),
data.begin() );
load_time = stop_clock();
solve_max();
max_time = stop_clock();
solve_min();
min_time = stop_clock();

cout << "Execution time: " << endl;
cout << "Data load " << load_time << endl;
cout << "Minimum " << min_time << endl;
cout << "Maximum " << max_time << endl;
cout << "====================" << endl;
cout << "Total " << load_time + min_time + max_time << endl;
}

[gugoXX]

ResMin : 75ms Dist: 0.000366913641159708 Indexes: 22503-72620
ResMax : 6ms Dist: 1.70445415447397 Indexes: 72832-48483
81ms


Eccomi.
Il minimo e' esatto, il massimo e' greedy.
Nella soluzione di ieri anche il massimo era esatto, ma i tempi oltre 2 secondi.
Il power della greedy si fa sentire :)
Pero' appunto e' greedy. Non c'e' garanzia che la soluzione sia quella giusta, anche se un risultato molto buono dovrebbe comparire il piu' delle volte.
Ora organizzo il codice e posto.

Qui il, file da 1milione di punti

ResMin : 2175ms Dist: 0.000353266507043159 Indexes: 830311-849536
ResMax : 43ms Dist: 1.72516305687739 Indexes: 575196-111795
2218ms


PS:Ho visto ora marco.r
La ricerca del minimo penso sia identica alla tua, e quella del massimo molto simile.
In quella del minimo hai pero' commesso un errore. Cosa capita se i 2 punti piu' vicini sono proprio a cavallo di uno dei miniscatoli?

[gugoXX]

Ecco


static void Main(string[] args)
{
Universe dc = new Universe();
dc.Load(@"C:\temp\Coords.dat");
Stopwatch sw = new Stopwatch();
sw.Start();
ResearchMin clus = new ResearchMin(dc.Allpoints, dc.LatoCubo);
var resmin = clus.Research();
long mintime = sw.ElapsedMilliseconds;
ResearchMax2 rmax = new ResearchMax2(dc.Allpoints, dc.LatoCubo);
var resmax = rmax.Research();

sw.Stop();
Console.WriteLine("ResMin : {0}ms {1}", mintime, resmin);
Console.WriteLine("ResMax : {0}ms {1}", sw.ElapsedMilliseconds-mintime, resmax);
Console.WriteLine("{0}ms",sw.ElapsedMilliseconds);
Console.ReadKey();
}


public class Point3d
{
public double x;
public double y;
public double z;
public Point3d(double ix,double iy,double iz)
{
x = ix;
y = iy;
z = iz;
}

public static Random rnd = new Random(18);
public static Point3d Rand(double dim)
{
return new Point3d(rnd.NextDouble() * dim, rnd.NextDouble() * dim, rnd.NextDouble() * dim);
}

public static double operator -(Point3d d1, Point3d d2)
{
double dx=d2.x-d1.x;
double dy=d2.y-d1.y;
double dz=d2.z-d1.z;
return (dx * dx + dy * dy + dz * dz);
}

public override string ToString()
{
return string.Format("{0} {1} {2}", x, y, z);
}
}


public class Result
{
public Result(double init)
{
dis = init;
}

public Result(double p_dis,int i1,int i2)
{
dis = p_dis;
index1 = i1;
index2 = i2;
}

public double dis;
public int index1, index2;

public override string ToString()
{
return string.Format("Dist: {0} Indexes: {1}-{2}", Math.Sqrt(dis), index1, index2);
}

public static Result Greatest(Result r1, Result r2)
{
if (r1.dis > r2.dis) return r1;
return r2;
}
}

public class Point3dIndex : Point3d
{
public Point3dIndex(double x, double y, double z,int Idx)
:base(x,y,z)
{
Index = Idx;
}

public int Index;
}

public class ResearchMin
{
List<Point3d> allpoints;
int SplitFactor;
double LatoCubo;

public ResearchMin(List<Point3d> input, double LCubo)
{
allpoints = input;
LatoCubo = LCubo;
}

List<Point3dIndex>[, ,] Zoner;

public Result Research()
{
int ln = allpoints.Count;
if (ln == 1){
return new Result(0,0,0);
}
if (ln == 2)
{
double dis = allpoints[1] - allpoints[2];
Result rr = new Result(dis,0,1);
}

Result res = new Result(double.MaxValue);
double SplitFactorD = (int)(2 * Math.Sqrt(Math.Sqrt(allpoints.Count)));
int SplitFactor = (int)SplitFactorD;

RsMin(SplitFactor, res, false);
RsMin(SplitFactor, res, true);

if (res.dis > (LatoCubo / SplitFactorD))
throw new Exception("Minicubi troppo piccoli (Pochi punti, soluzione classica)");

return res;
}

private Result RsMin(int sp, Result res,bool Sfasa)
{
SplitFactor = sp;

double DimCella = LatoCubo/(double)SplitFactor;

// Initialize Cluster
int SpFSf = SplitFactor + (Sfasa?1:0);
Zoner = new List<Point3dIndex>[SpFSf, SpFSf, SpFSf];
for (int x = 0; x < SpFSf; x++)
{
for (int y = 0; y < SpFSf; y++)
{
for (int z = 0; z < SpFSf; z++)
{
Zoner[x, y, z] = new List<Point3dIndex>();
}
}
}

double sfasa=(Sfasa) ? (DimCella/2) : 0;
// Clusterize
int ln = allpoints.Count;
for (int t=0;t<ln;t++)
{
Point3d p3d = allpoints[t];
Point3dIndex p3di = new Point3dIndex(p3d.x, p3d.y, p3d.z, t);

int nx = Normalize(p3d.x+sfasa, LatoCubo, SplitFactor);
int ny = Normalize(p3d.y+sfasa, LatoCubo, SplitFactor);
int nz = Normalize(p3d.z+sfasa, LatoCubo, SplitFactor);

Zoner[nx, ny, nz].Add(p3di);
}

for (int x = 0; x < SpFSf; x++)
{
for (int y = 0; y < SpFSf; y++)
{
for (int z = 0; z < SpFSf; z++)
{
MinWithinList1(Zoner[x, y, z], res);
}
}
}

return res;
}

private void MinWithinList1(List<Point3dIndex> src,Result CurrentResult)
{
int ln = src.Count;
for (int t = 0; t < ln-1; t++)
{
Point3dIndex p3di1=src[t];
for (int u = t+1; u < ln; u++)
{
Point3dIndex p3di2 = src[u];
double dist = p3di2 - p3di1;

if (dist < CurrentResult.dis)
{
CurrentResult.dis = dist;
CurrentResult.index1 = p3di1.Index;
CurrentResult.index2 = p3di2.Index;
}
}
}
}

private int Normalize(double val, double max, int SplitFactor)
{
return (int)((val / max) * (double)SplitFactor);
}
}

public class ResearchMax2
{
private class rm2
{
public rm2(Point3d rif)
{
riferimento = rif;
}
public Point3dIndex punto;
public Point3d riferimento;
public double cdist = double.MaxValue;
}

List<Point3d> allpoints;
double LatoCubo;

public ResearchMax2(List<Point3d> ap, double lc)
{
allpoints = ap;
LatoCubo = lc;
}

rm2 cmmm = new rm2(new Point3d(0, 0, 0));
rm2 cmmp = new rm2(new Point3d(0, 0, 1));
rm2 cmpm = new rm2(new Point3d(0, 1, 0));
rm2 cmpp = new rm2(new Point3d(0, 1, 1));
rm2 cpmm = new rm2(new Point3d(1, 0, 0));
rm2 cpmp = new rm2(new Point3d(1, 0, 1));
rm2 cppm = new rm2(new Point3d(1, 1, 0));
rm2 cppp = new rm2(new Point3d(1, 1, 1));

public Result Research()
{
int ln = allpoints.Count;
for (int i = 0; i < ln; i++)
{
Point3d p3d = allpoints[i];
double x = p3d.x;
double y = p3d.y;
double z = p3d.z;
double LMez = LatoCubo / 2;
rm2 pref = null;
if (x < LMez && y < LMez && z < LMez) pref = cmmm;
else if (x < LMez && y < LMez && z >= LMez) pref = cmmp;
else if (x < LMez && y >= LMez && z < LMez) pref = cmpm;
else if (x < LMez && y >= LMez && z >= LMez) pref = cmpp;
else if (x >= LMez && y < LMez && z < LMez) pref = cpmm;
else if (x >= LMez && y < LMez && z >= LMez) pref = cpmp;
else if (x >= LMez && y >= LMez && z < LMez) pref = cppm;
else if (x >= LMez && y >= LMez && z >= LMez) pref = cppp;

double dst = p3d - pref.riferimento;
if (dst < pref.cdist)
{
pref.cdist = dst;
pref.punto = new Point3dIndex(p3d.x, p3d.y, p3d.z, i);
}
}

double d1 = cmmm.punto - cppp.punto;
Result r1 = new Result(d1, cmmm.punto.Index, cppp.punto.Index);

double d2 = cmmp.punto - cppm.punto;
Result r2 = new Result(d2, cmmp.punto.Index, cppm.punto.Index);

double d3 = cmpm.punto - cpmp.punto;
Result r3 = new Result(d3, cmpm.punto.Index, cpmp.punto.Index);

double d4 = cmpp.punto - cpmm.punto;
Result r4 = new Result(d4, cmpp.punto.Index, cpmm.punto.Index);

Result CurrentResult = Result.Greatest(r1, r2);
CurrentResult = Result.Greatest(CurrentResult, r3);
CurrentResult = Result.Greatest(CurrentResult, r4);

return CurrentResult;
}
}

[marco.r]

PS:Ho visto ora marco.r
La ricerca del minimo penso sia identica alla tua, e quella del massimo molto simile.
In quella del minimo hai pero' commesso un errore. Cosa capita se i 2 punti piu' vicini sono proprio a cavallo di uno dei miniscatoli?
Confronto sempre ogni miniscatolo (!) (:D) con i sette adiacenti successivi per evitare questo. Ripeto un po' di conti ma mi evita l'errore.
L'algoritmo non e' ottimale, a naso quello di repne scala meglio, pero' il mio e' talmente "lineare" che fino a un milione di elementi e' ancora piu' veloce.

[gugoXX]

Confronto sempre ogni miniscatolo (!) (:D) con i sette adiacenti successivi per evitare questo. Ripeto un po' di conti ma mi evita l'errore.

Giusto. Ma i miniscatoli 3d vicini sono 26 (3^3-1), o forse non ho capito bene.
Comunque la mia e' simile. Faccio una prima passata con i miniscatoli come i tuoi, poi sfaso i miniscatoli un po' e ripasso, e poi un altro po' e ripasso di nuovo. In pratica faccio in modo che qualsiasi eventuale minimo accavallato sia compreso completamente in uno scatolo almeno una volta in una delle 3 configurazioni.

[marco.r]

Ecco i risultati sulla mia macchina
(portatile core 2 duo, con flag -march=core2 -msse3 -O3)

Maximum distance is 1.72516 between [575196,111795]
Minimum distance is 0.000353267 between [830311,849536]
Execution time:
Data load 2.16
Minimum 1.8
Maximum 0.03
====================
Total 3.99
mrighele@sasuke 8 % ./repne Coords.dat
Distanza minima: 0.000353266507043 - [830311][849536]
Distanza massima: 1.725163056877385 - [575196][111795]

*** Tempo per lettura dati: 1.670 secondi

*** Tempo per xxxxxxxxxxxxx: 3.020 secondi
*** Tempo per ricerca minimo: 1.900 secondi
*** Tempo per ricerca massimo: 0.640 secondi

*** Tempo totale soluzione: 5.560 secondi

La cosa carina e' che il mio e' piu' veloce nel minimo, quello di repne nel massimo.

[marco.r]

Giusto. Ma i miniscatoli 3d vicini sono 26 (3^3-1), o forse non ho capito bene.
Comunque la mia e' simile. Faccio una prima passata con i miniscatoli come i tuoi, poi sfaso i miniscatoli un po' e ripasso, e poi un altro po' e ripasso di nuovo. In pratica faccio in modo che qualsiasi eventuale minimo accavallato sia compreso completamente in uno scatolo almeno una volta in una delle 3 configurazioni.
Hum, hai in parte ragione, nel senso salto alcune coppie. In realta' non mi servono tutte perche' alcune le verifico quando controllo altri cubi, comunque ne devo controllare 19 :S. Il tempo presumibilmente aumentera' di un 150% :stordita:

[marco.r]

Giusto. Ma i miniscatoli 3d vicini sono 26 (3^3-1), o forse non ho capito bene.
Comunque la mia e' simile. Faccio una prima passata con i miniscatoli come i tuoi, poi sfaso i miniscatoli un po' e ripasso, e poi un altro po' e ripasso di nuovo. In pratica faccio in modo che qualsiasi eventuale minimo accavallato sia compreso completamente in uno scatolo almeno una volta in una delle 3 configurazioni.

Anche, buona idea

[rеpne scasb]

■

[Vincenzo1968]

Mbè?

se pò avè la spiegazzione de l'argoritmo? Io ho trovato quarche cosa sul web e drento quarche libro, ma gli esempi tratteno tutti er caso bidimensionale.
Gradirei una quarche spiegazzione, co' l'ausilio de quarche grafico, de quarche disegnino, in 3D.

Grazzie :)

:bimbo:

[gugoXX]

Mbè?

se pò avè la spiegazzione de l'argoritmo? Io ho trovato quarche cosa sul web e drento quarche libro, ma gli esempi tratteno tutti er caso bidimensionale.
Gradirei una quarche spiegazzione, co' l'ausilio de quarche grafico, de quarche disegnino, in 3D.

Grazzie :)

:bimbo:

Beh, l'algoritmo per trovare il minimo l'avevo gia' postato, e lo riporto:
Suddivido il cubo in miniscatoli, e cerco il minimo dentro ciascuno dei miniscatoli.
Essendo che sto cercando il minimo infatti non mi serve andare a calcolare le distanze relative tra i punti di un miniscatolo e quelli di un altro miniscatolo, in quanto assumo che il minimo che trovero' sara' piu' piccolo del lato del miniscatolo. Se cosi' non fosse dovrei aumentare il lato di ciascun miniscatolo.
C'e' il problema dell'accavallamento: se i 2 punti piu' vicini fossero proprio a cavallo di due miniscatoli, rischierei di perderli.
Quindi faccio una prima passata con i miniscatoli come detto, poi sfaso i miniscatoli un po' e ripasso, e poi un altro po' e ripasso di nuovo. In pratica faccio in modo che qualsiasi eventuale minimo accavallato sia compreso completamente in uno scatolo almeno una volta in una delle 3 configurazioni.


L'algoritmo greedy per il massimo e' ovviamente semplice (come giusto che sia per un greedy).
Per ciascuno degli 8 vertici del cubo cerco il punto piu' vicino.
Calcolo poi la distanza tra il punto cosi' trovato di un vertice e quello relativo al vertice opposto. Tengo la coppia di punti piu' distante.
Ovviamente ci sono casi in cui non funziona, ma e' mediamente corretto, sbaglia quando i punti sono pochi o quando sono distribuiti in modo "strano"

[Vincenzo1968]

Ahò Gugo,

Grazzie assai ;)
Gradirei però, anche una sorta de pseudo-codice e, soprattutto, quarche bel disegnino(come ner caso del a spiegazione che hai postato ner contest sul sotto-array de somma massima).

;)

[Vincenzo1968]

Ahò,

ce so' riuscito(almeno in parte), li mortan guerieri!

file da 100.000 punti
http://www.guidealgoritmi.it/images/ImgForums/Contest08Soluz01OptTemp2.jpg

file da 1.000.000 de punti
http://www.guidealgoritmi.it/images/ImgForums/Contest08Soluz01OptTemp3.jpg


La ricerca der minimo è coretta. Devo aggiustà er tiro pe la ricerca der massimo.
So' forti st'alberi tridimensionali ;)

:bimbo:

P.S Er codice lo posto quanno che ho risorto anco er problema der massimo.

[Vincenzo1968]

Ahò,

riesco a batte' Repne sui tempi pe la ricerca der minimo.

Questi so' i tempi de repne sulla mia macchina(file da 1.000.000 di punti):
http://www.guidealgoritmi.it/images/ImgForums/Contest08ASoluzRepne2.jpg

Mò vojo provà a prenne i tempi della soluzione de Gugo.

:bimbo:

[Vincenzo1968]

Er programma de Marco.r me va in crash sul file grosso(quelo da 1.000.000 de punti).
Er programma de Gugo, invece, nun se compila perchè nun riconosce er tipo 'Universe'.

Posto er codice mio anche se sbaglia er calcolo de la distanza massima:


#include <stdio.h>
#include <math.h>
#include <time.h>
#include <malloc.h>

#define BUFFER_SIZE 4096
#define MAX_STACK 100

typedef double ElementType;
typedef ElementType ItemType[3];

typedef struct tagPunti
{
ItemType data;
int index;
} Punti;

typedef struct tagRes
{
Punti P1;
Punti P2;
double Distanza;
int index1;
int index2;
} Res;

typedef struct tagArray
{
int dimensione;
double dmax;
Punti *m;
} Array;

typedef struct tagTree
{
Punti P;
//struct tagTree *father;
struct tagTree *left;
struct tagTree *right;
} Tree;

Tree *NewNode(ItemType key, int index);
Tree *InsertNode(Tree *node, ItemType key, int index);
void FreeTree(Tree *head);


typedef enum tagStati
{
S_ERROR = -1, S0 = 0, S1, S2, S3, S4, S5, S6, S7, S8, S9
} Stati;


// La macro seguente restituisce il numero degli elementi di un array.
// Per esempio, data la seguente dichiarazione:
//
// int array[] = {8, 5, 13, 2, 1};
//
// ARRAY_SIZE(array) restituisce 5

#define ARRAY_SIZE(Array) (sizeof(Array) / sizeof((Array)[0]))

Stati DFA(const char *szFileName, Array *pArray)
{
Stati stato = S0;
FILE *fp;
unsigned char buffer[BUFFER_SIZE];
int numblocks;
int numread;
unsigned char c;
int k, x, j;
int riga, colonna;
char szNum[256];
double num;

unsigned char byteCR = 0xD; // Carriage Return
unsigned char byteLF = 0xA; // Line Feed

fp = fopen(szFileName, "rb");
if ( fp == NULL )
{
printf("Errore nell'apertura del file %s\n", szFileName);
return S_ERROR;
}

if ( fseek(fp, 0, SEEK_END) )
return 0;

numblocks = ftell(fp)/BUFFER_SIZE;
if ( numblocks == 0 )
{
numblocks = 1;
}
else
{
if ( ftell(fp) % BUFFER_SIZE != 0 )
numblocks++;
}

fseek(fp, 0, SEEK_SET);
numread = fread(buffer, 1, BUFFER_SIZE, fp);
if (numread == 0)
{
fclose(fp);
printf("Errore 1 nella lettura del file %s\n", szFileName);
return S_ERROR;
}

pArray->dimensione = 0;
k = 0;
c = *(buffer + k++);
while ( 1 )
{
if ( c >= '0' && c <= '9' )
{
pArray->dimensione = pArray->dimensione * 10 + c - '0';
}
else if ( c == ' ' || c == '\t' )
{
j = 0;
while ( c != byteLF )
{
c = *(buffer + k++);
if ( c >= '0' && c <= '9' )
szNum[j++] = c;
}
szNum[j] = '\0';
pArray->dmax = atof(szNum);
break;
}
else if ( c == '\n' )
{
break;
}

c = *(buffer + k++);
}

pArray->m = (Punti*)malloc(pArray->dimensione*sizeof(Punti));
if ( !pArray->m )
{
printf("Memoria insufficiente.\n");
fclose(fp);
return S_ERROR;
}

riga = colonna = 0;
num = 0;
x = j = 0;
while ( x < numblocks )
{
c = *(buffer + k++);

if ( c == byteLF )
{
pArray->m[riga].index = riga;
riga++;
colonna = 0;
}

switch (stato)
{
case S0:
j = 0;
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S2;
}
else if ( c == '.' )
{
szNum[j++] = c;
stato = S3;
}
else if ( c == '-' || c == '+' )
{
szNum[j++] = c;
stato = S1;
}
else if ( c == ' ' || c == '\t' )
{
while ( c == ' ' || c == '\t' )
{
c = *(buffer + k++);
if ( k >= numread )
break;
}
k--;
}
else if ( c == byteCR || c == byteLF )
{
// nada
}
else
{
printf("Errore: stato S0; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S1:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S2;
}
else if ( c == '.' )
{
szNum[j++] = c;
stato = S3;
}
else
{
printf("Errore: stato S1; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S2:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
}
else if ( c == '.' )
{
szNum[j++] = c;
stato = S4;
}
else if ( c == 'E' || c == 'e' )
{
szNum[j++] = c;
stato = S5;
}
else if ( c == ' ' || c == '\t' || c == byteCR || c == byteLF )
{
szNum[j] = '\0';
pArray->m[riga].data[colonna++] = atof(szNum);
stato = S0;
}
else
{
printf("Errore: stato S2; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S3:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S6;
}
else
{
printf("Errore: stato S3; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S4:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S6;
}
else if ( c == 'E' || c == 'e' )
{
szNum[j++] = c;
stato = S7;
}
else if ( c == ' ' || c == '\t' || c == byteCR || c == byteLF )
{
szNum[j] = '\0';
pArray->m[riga].data[colonna++] = atof(szNum);
stato = S0;
}
else
{
printf("Errore: stato S4; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S5:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S9;
}
else if ( c == '-' || c == '+' )
{
szNum[j++] = c;
stato = S8;
}
else
{
printf("Errore: stato S5; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S6:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S6;
}
else if ( c == 'E' || c == 'e' )
{
szNum[j++] = c;
stato = S7;
}
else if ( c == ' ' || c == '\t' || c == byteCR || c == byteLF )
{
szNum[j] = '\0';
pArray->m[riga].data[colonna++] = atof(szNum);
stato = S0;
}
else
{
printf("Errore: stato S6; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S7:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S9;
}
else if ( c == '-' || c == '+' )
{
szNum[j++] = c;
stato = S8;
}
else
{
printf("Errore: stato S7; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S8:
case S9:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S9;
}
else if ( c == ' ' || c == '\t' || c == byteCR || c == byteLF )
{
szNum[j] = '\0';
pArray->m[riga].data[colonna++] = atof(szNum);
stato = S0;
}
else
{
printf("Errore: stato S9; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;
}

if ( k >= numread )
{
numread = fread(buffer, 1, BUFFER_SIZE, fp);
k = 0;
x++;
}
}

fclose(fp);

return stato;
}

Tree *NewNode(ItemType key, int index)
{
Tree *r;

r = (Tree *) malloc(sizeof(Tree));

if(!r)
{
printf("Memoria insufficiente.\n");
return NULL;
}

r->P.data[0] = key[0];
r->P.data[1] = key[1];
r->P.data[2] = key[2];
r->P.index = index;
//r->father = NULL;
r->left = NULL;
r->right = NULL;

return r;
}

Tree *InsertNode(Tree *node, ItemType key, int index)
{
Tree *pRadice = NULL;
int Level = 0;

if( !node )
{
node = NewNode(key, index);
return node;
}

pRadice = node;

while( 1 )
{
if ( key[Level] < node->P.data[Level] )
{
if ( !node->left )
{
node->left = NewNode(key, index);
//node->left->father = node;
break;
}
node = node->left;
}
else
{
if ( !node->right )
{
node->right = NewNode(key, index);
//node->right->father = node;
break;
}
node = node->right;
}
Level++;
if ( Level > 2 )
Level = 0;
}

node = pRadice;

return node;
}

void FreeTree(Tree *head)
{
Tree *temp1, *temp2;

Tree *stack[MAX_STACK];
int top;

top = 0;

if ( !head )
return;

temp1 = temp2 = head;

while ( temp1 != NULL )
{
for(; temp1->left != NULL; temp1 = temp1->left)
stack[top++] = temp1;

while ( (temp1 != NULL) && (temp1->right == NULL || temp1->right == temp2) )
{
temp2 = temp1;
free(temp2);
if ( top == 0 )
return;
temp1 = stack[--top];
}

stack[top++] = temp1;
temp1 = temp1->right;
}
}

void TreeToArray(Tree *head, Array *pArray)
{
Tree *temp;

Tree *stack[MAX_STACK];
int top = 0;
int k = 0;

if ( !head )
return;

temp = head;

stack[top++] = temp; // Push

while ( top != 0 )
{
temp = stack[--top]; // Pop

//printf("%d -> (%lf,%lf,%lf)\n", temp->index, temp->data[0], temp->data[1], temp->data[2]);
pArray->m[k++] = temp->P;

if ( temp->right != NULL )
stack[top++] = temp->right;
if ( temp->left != NULL )
stack[top++] = temp->left;
}
}

void Calcola(Punti P[], int dim, Res *pMin, Res *pMax)
{
int x;
double distanza, distanzaMin, distanzaMax;
double diffA, diffB, diffC;

distanzaMin = distanzaMax = -1;

for ( x = 1; x < dim; x++ )
{

diffA = P[x-1].data[0] - P[x].data[0];
diffB = P[x-1].data[1] - P[x].data[1];
diffC = P[x-1].data[2] - P[x].data[2];
distanza = diffA*diffA + diffB*diffB + diffC*diffC;
if ( distanza < distanzaMin || distanzaMin == -1 )
{
distanzaMin = distanza;
pMin->P1 = P[x-1];
pMin->P2 = P[x];
pMin->index1 = P[x-1].index;
pMin->index2 = P[x].index;
}
/*
if ( distanza > distanzaMax )
{
distanzaMax = distanza;
pMax->P1 = P[x-1];
pMax->P2 = P[x];
pMax->index1 = P[x-1].index;
pMax->index2 = P[x].index;
}
*/
}

diffA = P[0].data[0] - P[dim-1].data[0];
diffB = P[0].data[1] - P[dim-1].data[1];
diffC = P[0].data[2] - P[dim-1].data[2];
distanzaMax = diffA*diffA + diffB*diffB + diffC*diffC;
pMax->P1 = P[0];
pMax->P2 = P[dim-1];
pMax->index1 = P[0].index;
pMax->index2 = P[dim-1].index;

pMin->Distanza = sqrt( distanzaMin );
pMax->Distanza = sqrt( distanzaMax );
}

int main()
{
Stati stato;
Array myArray;
Res r1, r2;
Tree *pTree = NULL;
int k;

clock_t c_start, c_end;

char *szFileName = "C:\\Scaricamenti\\Temp\\Gugo\\Contest 08 - Compattamento Barionico\\Coords2\\Coords.dat";

c_start = clock();
stato = DFA(szFileName, &myArray);
c_end = clock();
printf("\nTempo impiegato per caricamento dati -> %5.5f secondi\n", (double)(c_end - c_start) / CLOCKS_PER_SEC);

c_start = clock();
for(k = 0; k < myArray.dimensione; k++)
pTree = InsertNode(pTree, myArray.m[k].data, myArray.m[k].index);
TreeToArray(pTree, &myArray);
Calcola(myArray.m, myArray.dimensione, &r1, &r2);
c_end = clock();

printf("\nDistanza Min: P[%d] - P[%d] -> %.15lf\n", r1.index1, r1.index2, r1.Distanza);
printf("\nDistanza Max: P[%d] - P[%d] -> %.15lf\n", r2.index1, r2.index2, r2.Distanza);
printf("\nTempo impiegato per la ricerca -> %5.5f secondi\n", (double)(c_end - c_start) / CLOCKS_PER_SEC);

free(myArray.m);
FreeTree(pTree);

return 0;
}


:bimbo:

[banryu79]

@Vincenzo1968: ma per quale stramaledetto motivo scrivi come se parlassi in romanaccio, nun se capisce na mazza :D

[Vincenzo1968]

@Vincenzo1968: ma per quale stramaledetto motivo scrivi come se parlassi in romanaccio, nun se capisce na mazza :D

No savaria dirghe :boh:
Forse perché, in questo periodo, sto rileggendo il Pasticciaccio di Gadda(El me ga ipnotisà).

:bimbo:

x Gugo: ohé, manca il pezzo di codice relativo alla classe 'Universe'. Quanno te decidi a postallo? :read:

:bimbo:

[^TiGeRShArK^]

No savaria dirghe :boh:
Forse perché, in questo periodo, sto rileggendo il Pasticciaccio di Gadda(El me ga ipnotisà).

:bimbo:

x Gugo: ohé, manca il pezzo di codice relativo alla classe 'Universe'. Quanno te decidi a postallo? :read:

:bimbo:

ahò, non ce stà 'r codice de l'universo. Ce lo voi postà? :O

:asd:

[Vincenzo1968]

ahò, non ce stà 'r codice de l'universo. Ce lo voi postà? :O

:asd:

:D

Vedi com'è divertente? Il dialetto romanesco non è poi così difficile da capire(grazie anche all'uso che se ne fa in letteratura e, soprattutto, nel cinema) ;)

[^TiGeRShArK^]

:D

Vedi com'è divertente? Il dialetto romanesco non è poi così difficile da capire(grazie anche all'uso che se ne fa in letteratura e, soprattutto, nel cinema) ;)

io più che altro sono stato poco + di un anno a lavorare a roma e con i miei colleghi ormai m'ero abituato :asd:

[gugoXX]

Ecco, ecco...


public class Universe
{
public List<Point3d> Allpoints;

public Universe()
{
Allpoints = new List<Point3d>();
}

public void Load(string fname)
{
StreamReader sr = new StreamReader(fname);
string stt = sr.ReadLine();
string[] first = stt.Split(new char[] { ' ' }, StringSplitOptions.RemoveEmptyEntries);
int howmany = int.Parse(first[0]);
LatoCubo = double.Parse(first[1]);
string line;
while ((line = sr.ReadLine()) != null)
{
string[] coords = line.Split(new char[] { ' ' }, StringSplitOptions.RemoveEmptyEntries);
Allpoints.Add(new Point3d(double.Parse(coords[0]),double.Parse(coords[1]),double.Parse(coords[2])));
}
sr.Close();
}

public void Save(string fname)
{
StreamWriter sw = new StreamWriter(fname);
sw.WriteLine("{0} {1}", Allpoints.Count, LatoCubo);
foreach (Point3d p3d in Allpoints)
{
sw.WriteLine(p3d.ToString());
}
sw.Close();
}

public double LatoCubo;
public void Generate(int npoints,double dim)
{
LatoCubo = dim;
Allpoints.Clear();
for (int t = 0; t < npoints; t++)
{
Allpoints.Add(Point3d.Rand(dim));
}
}
}

[Vincenzo1968]

Ecco, ecco...


public class Universe
{
public List<Point3d> Allpoints;

public Universe()
{
Allpoints = new List<Point3d>();
}

public void Load(string fname)
{
StreamReader sr = new StreamReader(fname);
string stt = sr.ReadLine();
string[] first = stt.Split(new char[] { ' ' }, StringSplitOptions.RemoveEmptyEntries);
int howmany = int.Parse(first[0]);
LatoCubo = double.Parse(first[1]);
string line;
while ((line = sr.ReadLine()) != null)
{
string[] coords = line.Split(new char[] { ' ' }, StringSplitOptions.RemoveEmptyEntries);
Allpoints.Add(new Point3d(double.Parse(coords[0]),double.Parse(coords[1]),double.Parse(coords[2])));
}
sr.Close();
}

public void Save(string fname)
{
StreamWriter sw = new StreamWriter(fname);
sw.WriteLine("{0} {1}", Allpoints.Count, LatoCubo);
foreach (Point3d p3d in Allpoints)
{
sw.WriteLine(p3d.ToString());
}
sw.Close();
}

public double LatoCubo;
public void Generate(int npoints,double dim)
{
LatoCubo = dim;
Allpoints.Clear();
for (int t = 0; t < npoints; t++)
{
Allpoints.Add(Point3d.Rand(dim));
}
}
}


Me dà erore: indice oltre i limiti della matrice. :mad:

[gugoXX]

Me dà erore: indice oltre i limiti della matrice. :mad:

Dove? Il file di input e' secondo specifiche?

[Vincenzo1968]

Dove? Il file di input e' secondo specifiche?

Ho provato sia col file tuo(quello da 100.000 punti )che con quello di repne(file da 1.000.000 di punti).
A runtime ho quell'errore. Te posso posta' er screenshotte.

:bimbo:

[Vincenzo1968]

Eqque qua:

http://www.guidealgoritmi.it/images/ImgForums/Contest08ErroreGugo.jpg

[Vincenzo1968]

In omaggio al dialetto romanesco, ecco una poesia di Trilussa(e ringrazio Gugo per avermela segnalata):


LA POESIA DELL'UCELLETTO

Era d’Agosto e un povero uccelletto
ferito dallo fionda di un maschietto
andò per riposare l’ala offesa,
sulla finestra aperta di una chiesa.
Dalle tendine del confessionale
il parroco intravide l’animale
ma, pressato dal ministero urgente,
rimase intento a confessar la gente.
Mentre in ginocchio alcuni altri a sedere
dicevano i fedeli le preghiere,
una donna, notato l’uccelletto,
lo prese al caldo e se lo mise al petto.
D'un tratto un cinguettio ruppe il silenzio
e il prete a quel rumore
il ruolo abbandonò di confessore.
Scuro in viso peggio della pece,
s'arrampicò sul pulpito e poi fece:
“Fratelli! Chi ha l’uccello per favore
esca fuori dal tempio del Signore!”
I maschi, un po’ stupiti a tal parole,
lenti s'accinsero ad alzar le suole,
ma il prete a quell'errore madornale
“Fermi” gridò “mi sono espresso male!
Rientrate tutti e statemi a sentire,
solo chi ha preso l’uccello deve uscire!”
A testa bassa, la corona in mano,
cento donne s'alzarono pian piano.
Ma mentre se ne andavano ecco allora che il parroco strillò:
“Sbagliate ancora, rientrate tutte quante figlie amate
che in non volevo dir quel che pensate!
Ecco, quello che ho detto torno a dire,
solo chi ha preso l'uccello deve uscire,
ma, mi rivolgo, non ci sia sorpresa,
soltanto a chi l'uccello l’ha preso in chiesa!”
Finì la frase e nello stesso istante
le monache s'alzarono tutte quante
e con il volto pieno di rossore
lasciavano la casa del Signore.
“O Santa Vergine!” esclamo il buon prete
“Fatemi la grazia se potete.
Poi senza fare rumore dico, piano piano
s'alzi soltanto chi ha l’uccello in mano!”
Una ragazza che col fidanzato s'era messa in un angolo appartato
sommessa mormorò con viso smorto
“Che ti dicevo, hai visto? Se n’è accorto!”


Qui è recitata da Andrea Bocelli:

http://it.youtube.com/watch?v=sYB7Ev5yLR0

:bimbo:

[Vincenzo1968]

Ho risolto il problema anche per la ricerca della distanza massima. Adesso il risultato è corretto ma l'algoritmo è lento; impiega più di trenta secondi sul file grosso(quello postato da repne, da un milone di record).

Per la ricerca della distanza minima, invece, i tempi risultano, incredibilmente, più veloci di quelli di repne :confused:

File da 100.000 records:
http://www.guidealgoritmi.it/images/ImgForums/Contest08SoluzFilePiccolo.jpg

File da 1.000.000 di records:
http://www.guidealgoritmi.it/images/ImgForums/Contest08SoluzFileGrosso.jpg

La mia macchina:

AMD Athlon(tm) 64 X2
Dual Core Processor 4800+
2.50 GHz
896 MB di RAM

Microsoft Windows XP Professional (32 bit)
Service Pack 3



#include <stdio.h>
#include <math.h>
#include <time.h>
#include <malloc.h>

#define BUFFER_SIZE 4096
#define MAX_STACK 100

typedef double ElementType;
typedef ElementType ItemType[3];

typedef struct tagPunti
{
ItemType data;
int index;
} Punti;

typedef struct tagRes
{
Punti P1;
Punti P2;
double Distanza;
int index1;
int index2;
} Res;

typedef struct tagArray
{
int dimensione;
double dmax;
Punti *m;
} Array;

typedef struct tagTree
{
Punti P;
//struct tagTree *father;
struct tagTree *left;
struct tagTree *right;
} Tree;

Tree *NewNode(ItemType key, int index);
Tree *InsertNode(Tree *node, ItemType key, int index);
void FreeTree(Tree *head);


typedef enum tagStati
{
S_ERROR = -1, S0 = 0, S1, S2, S3, S4, S5, S6, S7, S8, S9
} Stati;


// La macro seguente restituisce il numero degli elementi di un array.
// Per esempio, data la seguente dichiarazione:
//
// int array[] = {8, 5, 13, 2, 1};
//
// ARRAY_SIZE(array) restituisce 5

#define ARRAY_SIZE(Array) (sizeof(Array) / sizeof((Array)[0]))

Stati DFA(const char *szFileName, Array *pArray)
{
Stati stato = S0;
FILE *fp;
unsigned char buffer[BUFFER_SIZE];
int numblocks;
int numread;
unsigned char c;
int k, x, j;
int riga, colonna;
char szNum[256];
double num;

unsigned char byteCR = 0xD; // Carriage Return
unsigned char byteLF = 0xA; // Line Feed

fp = fopen(szFileName, "rb");
if ( fp == NULL )
{
printf("Errore nell'apertura del file %s\n", szFileName);
return S_ERROR;
}

if ( fseek(fp, 0, SEEK_END) )
return 0;

numblocks = ftell(fp)/BUFFER_SIZE;
if ( numblocks == 0 )
{
numblocks = 1;
}
else
{
if ( ftell(fp) % BUFFER_SIZE != 0 )
numblocks++;
}

fseek(fp, 0, SEEK_SET);
numread = fread(buffer, 1, BUFFER_SIZE, fp);
if (numread == 0)
{
fclose(fp);
printf("Errore 1 nella lettura del file %s\n", szFileName);
return S_ERROR;
}

pArray->dimensione = 0;
k = 0;
c = *(buffer + k++);
while ( 1 )
{
if ( c >= '0' && c <= '9' )
{
pArray->dimensione = pArray->dimensione * 10 + c - '0';
}
else if ( c == ' ' || c == '\t' )
{
j = 0;
while ( c != byteLF )
{
c = *(buffer + k++);
if ( c >= '0' && c <= '9' )
szNum[j++] = c;
}
szNum[j] = '\0';
pArray->dmax = atof(szNum);
break;
}
else if ( c == '\n' )
{
break;
}

c = *(buffer + k++);
}

pArray->m = (Punti*)malloc(pArray->dimensione*sizeof(Punti));
if ( !pArray->m )
{
printf("Memoria insufficiente.\n");
fclose(fp);
return S_ERROR;
}

riga = colonna = 0;
num = 0;
x = j = 0;
while ( x < numblocks )
{
c = *(buffer + k++);

if ( c == byteLF )
{
pArray->m[riga].index = riga;
riga++;
colonna = 0;
}

switch (stato)
{
case S0:
j = 0;
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S2;
}
else if ( c == '.' )
{
szNum[j++] = c;
stato = S3;
}
else if ( c == '-' || c == '+' )
{
szNum[j++] = c;
stato = S1;
}
else if ( c == ' ' || c == '\t' )
{
while ( c == ' ' || c == '\t' )
{
c = *(buffer + k++);
if ( k >= numread )
break;
}
k--;
}
else if ( c == byteCR || c == byteLF )
{
// nada
}
else
{
printf("Errore: stato S0; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S1:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S2;
}
else if ( c == '.' )
{
szNum[j++] = c;
stato = S3;
}
else
{
printf("Errore: stato S1; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S2:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
}
else if ( c == '.' )
{
szNum[j++] = c;
stato = S4;
}
else if ( c == 'E' || c == 'e' )
{
szNum[j++] = c;
stato = S5;
}
else if ( c == ' ' || c == '\t' || c == byteCR || c == byteLF )
{
szNum[j] = '\0';
pArray->m[riga].data[colonna++] = atof(szNum);
stato = S0;
}
else
{
printf("Errore: stato S2; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S3:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S6;
}
else
{
printf("Errore: stato S3; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S4:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S6;
}
else if ( c == 'E' || c == 'e' )
{
szNum[j++] = c;
stato = S7;
}
else if ( c == ' ' || c == '\t' || c == byteCR || c == byteLF )
{
szNum[j] = '\0';
pArray->m[riga].data[colonna++] = atof(szNum);
stato = S0;
}
else
{
printf("Errore: stato S4; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S5:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S9;
}
else if ( c == '-' || c == '+' )
{
szNum[j++] = c;
stato = S8;
}
else
{
printf("Errore: stato S5; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S6:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S6;
}
else if ( c == 'E' || c == 'e' )
{
szNum[j++] = c;
stato = S7;
}
else if ( c == ' ' || c == '\t' || c == byteCR || c == byteLF )
{
szNum[j] = '\0';
pArray->m[riga].data[colonna++] = atof(szNum);
stato = S0;
}
else
{
printf("Errore: stato S6; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S7:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S9;
}
else if ( c == '-' || c == '+' )
{
szNum[j++] = c;
stato = S8;
}
else
{
printf("Errore: stato S7; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S8:
case S9:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S9;
}
else if ( c == ' ' || c == '\t' || c == byteCR || c == byteLF )
{
szNum[j] = '\0';
pArray->m[riga].data[colonna++] = atof(szNum);
stato = S0;
}
else
{
printf("Errore: stato S9; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;
}

if ( k >= numread )
{
numread = fread(buffer, 1, BUFFER_SIZE, fp);
k = 0;
x++;
}
}

fclose(fp);

return stato;
}

Tree *NewNode(ItemType key, int index)
{
Tree *r;

r = (Tree *) malloc(sizeof(Tree));

if(!r)
{
printf("Memoria insufficiente.\n");
return NULL;
}

r->P.data[0] = key[0];
r->P.data[1] = key[1];
r->P.data[2] = key[2];
r->P.index = index;
//r->father = NULL;
r->left = NULL;
r->right = NULL;

return r;
}

Tree *InsertNode(Tree *node, ItemType key, int index)
{
Tree *pRadice = NULL;
int Level = 0;

if( !node )
{
node = NewNode(key, index);
return node;
}

pRadice = node;

while( 1 )
{
if ( key[Level] < node->P.data[Level] )
{
if ( !node->left )
{
node->left = NewNode(key, index);
//node->left->father = node;
break;
}
node = node->left;
}
else
{
if ( !node->right )
{
node->right = NewNode(key, index);
//node->right->father = node;
break;
}
node = node->right;
}
Level++;
if ( Level > 2 )
Level = 0;
}

node = pRadice;

return node;
}

void FreeTree(Tree *head)
{
Tree *temp1, *temp2;

Tree *stack[MAX_STACK];
int top;

top = 0;

if ( !head )
return;

temp1 = temp2 = head;

while ( temp1 != NULL )
{
for(; temp1->left != NULL; temp1 = temp1->left)
stack[top++] = temp1;

while ( (temp1 != NULL) && (temp1->right == NULL || temp1->right == temp2) )
{
temp2 = temp1;
free(temp2);
if ( top == 0 )
return;
temp1 = stack[--top];
}

stack[top++] = temp1;
temp1 = temp1->right;
}
}

void TreeToArrayMin(Tree *head, Array *pArray)
{
Tree *temp;

Tree *stack[MAX_STACK];
int top = 0;
int k = 0;

if ( !head )
return;

temp = head;

stack[top++] = temp; // Push

while ( top != 0 )
{
temp = stack[--top]; // Pop

//printf("%d -> (%lf,%lf,%lf)\n", temp->index, temp->data[0], temp->data[1], temp->data[2]);
pArray->m[k++] = temp->P;

if ( temp->right != NULL )
stack[top++] = temp->right;
if ( temp->left != NULL )
stack[top++] = temp->left;
}
}

void MakeArrayMax(Array *pArray, Array *pArrayMax1, Array *pArrayMax2)
{
int k, j, x;
double decimo;
double distanza;
double diffA, diffB, diffC;
double range1, range2;

decimo = pArray->dmax / 15.0;
range1 = pArray->dmax - decimo;
range2 = pArray->dmax*2.0 - decimo;

pArrayMax1->m = (Punti*)malloc(pArray->dimensione/2*sizeof(Punti));
if ( !pArrayMax1->m )
{
printf("Memoria insufficiente.\n");
return;
}

pArrayMax2->m = (Punti*)malloc(pArray->dimensione/2*sizeof(Punti));
if ( !pArrayMax2->m )
{
printf("Memoria insufficiente.\n");
free(pArrayMax1->m);
return;
}

j = x = 0;
for ( k = 0; k < pArray->dimensione; k++ )
{
diffA = pArray->m[k].data[0];
diffB = pArray->m[k].data[1];
diffC = pArray->m[k].data[2];
distanza = diffA*diffA + diffB*diffB + diffC*diffC;

if ( (distanza > range1) && (distanza < pArray->dmax) )
{
pArrayMax1->m[j++] = pArray->m[k];
/*
if ( pArray->m[k].index == 48483 )
{
printf("Eccolo -> 48483\n");
}
*/
}
else if ( distanza > range2 )
{
pArrayMax2->m[x++] = pArray->m[k];
/*
if ( pArray->m[k].index == 72832 )
{
printf("Eccolo -> 72832\n");
}
*/
}

}

pArrayMax1->dimensione = j;
pArrayMax1->dmax = pArray->dmax;

pArrayMax2->dimensione = x;
pArrayMax2->dmax = pArray->dmax;
}

void CalcolaMin(Punti P[], int dim, Res *pMin)
{
int x;
double distanza, distanzaMin;
double diffA, diffB, diffC;

distanzaMin = -1;

for ( x = 1; x < dim; x++ )
{
diffA = P[x-1].data[0] - P[x].data[0];
diffB = P[x-1].data[1] - P[x].data[1];
diffC = P[x-1].data[2] - P[x].data[2];
distanza = diffA*diffA + diffB*diffB + diffC*diffC;
if ( distanza < distanzaMin || distanzaMin == -1 )
{
distanzaMin = distanza;
pMin->P1 = P[x-1];
pMin->P2 = P[x];
pMin->index1 = P[x-1].index;
pMin->index2 = P[x].index;
}
}

pMin->Distanza = sqrt( distanzaMin );
}

void CalcolaMax(Punti P1[], Punti P2[], int dim1, int dim2, Res *pMax)
{
int x, y;
double distanza, distanzaMax;
double diffA, diffB, diffC;

distanzaMax = -1;

x = y = 0;

while( x < dim1 )
{
//y = x + 1;
y = 0;
while( y < dim2 )
{
//distanza = sqrt( pow(P[x].A - P[y].A, 2.0) + pow(P[x].B - P[y].B, 2.0) + pow(P[x].C - P[y].C, 2.0) );
diffA = P1[x].data[0] - P2[y].data[0];
diffB = P1[x].data[1] - P2[y].data[1];
diffC = P1[x].data[2] - P2[y].data[2];
distanza = diffA*diffA + diffB*diffB + diffC*diffC;

if ( distanza > distanzaMax )
{
distanzaMax = distanza;
pMax->P1 = P1[x];
pMax->P2 = P2[y];
pMax->index1 = P1[x].index;
pMax->index2 = P2[y].index;
}
++y;
}
++x;
}

pMax->Distanza = sqrt( distanzaMax );
}

int main()
{
Stati stato;
Array myArray, myArrayMax1, myArrayMax2;
Res r1, r2;
Tree *pTree = NULL;
int k;
double t1, t2;

clock_t c_start, c_end;

//char *szFileName = "C:\\Temp\\Coords\\Coords.dat";
char *szFileName = "C:\\Temp\\Coords2\\Coords.dat";

c_start = clock();
stato = DFA(szFileName, &myArray);
c_end = clock();
printf("\nTempo impiegato per caricamento dati -> %5.5f secondi\n", (double)(c_end - c_start) / CLOCKS_PER_SEC);

c_start = clock();
MakeArrayMax(&myArray, &myArrayMax1, &myArrayMax2);
for(k = 0; k < myArray.dimensione; k++)
pTree = InsertNode(pTree, myArray.m[k].data, myArray.m[k].index);
TreeToArrayMin(pTree, &myArray);
CalcolaMin(myArray.m, myArray.dimensione, &r1);
c_end = clock();
t1 = (double)(c_end - c_start) / CLOCKS_PER_SEC;

c_start = clock();
CalcolaMax(myArrayMax1.m, myArrayMax2.m, myArrayMax1.dimensione, myArrayMax2.dimensione, &r2);
c_end = clock();
t2 = (double)(c_end - c_start) / CLOCKS_PER_SEC;

printf("\nDistanza Min: P[%d] - P[%d] -> %.15lf\n", r1.index1, r1.index2, r1.Distanza);
printf("\nDistanza Max: P[%d] - P[%d] -> %.15lf\n", r2.index1, r2.index2, r2.Distanza);
printf("\nTempo impiegato per la ricerca della distanza minima -> %5.5f secondi\n", t1);
printf("\nTempo impiegato per la ricerca della distanza massima -> %5.5f secondi\n", t2);

free(myArray.m);
free(myArrayMax1.m);
free(myArrayMax2.m);
FreeTree(pTree);

return 0;
}


Volevo postare una classifica, ma:

La soluzione di Gugo mi dà questi errori a runtime:
http://www.guidealgoritmi.it/images/ImgForums/Contest08ErroreGugo.jpg

La soluzione di Marco.r, invece, va in crash.

:bimbo:

[gugoXX]

Rimetto qui il codice, ma dovrebbe essere uguale a prima


using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.IO;
using System.Diagnostics;

namespace Project6b
{
public static class MainClass
{
static void Main(string[] args)
{
Universe dc = new Universe();
dc.Load(@"C:\temp\Coords.dat");
Stopwatch sw = new Stopwatch();
sw.Start();
ResearchMin clus = new ResearchMin(dc.Allpoints, dc.LatoCubo);
var resmin = clus.Research();
long mintime = sw.ElapsedMilliseconds;
ResearchMax2 rmax = new ResearchMax2(dc.Allpoints, dc.LatoCubo);
var resmax = rmax.Research();

sw.Stop();
Console.WriteLine("ResMin : {0}ms {1}", mintime, resmin);
Console.WriteLine("ResMax : {0}ms {1}", sw.ElapsedMilliseconds - mintime, resmax);
Console.WriteLine("{0}ms", sw.ElapsedMilliseconds);
Console.ReadKey();
}
}


public class Point3d
{
public double x;
public double y;
public double z;
public Point3d(double ix, double iy, double iz)
{
x = ix;
y = iy;
z = iz;
}

public static Random rnd = new Random(18);
public static Point3d Rand(double dim)
{
return new Point3d(rnd.NextDouble() * dim, rnd.NextDouble() * dim, rnd.NextDouble() * dim);
}

public static double operator -(Point3d d1, Point3d d2)
{
double dx = d2.x - d1.x;
double dy = d2.y - d1.y;
double dz = d2.z - d1.z;
return (dx * dx + dy * dy + dz * dz);
}

public override string ToString()
{
return string.Format("{0} {1} {2}", x, y, z);
}
}


public class Result
{
public Result(double init)
{
dis = init;
}

public Result(double p_dis, int i1, int i2)
{
dis = p_dis;
index1 = i1;
index2 = i2;
}

public double dis;
public int index1, index2;

public override string ToString()
{
return string.Format("Dist: {0} Indexes: {1}-{2}", Math.Sqrt(dis), index1, index2);
}

public static Result Greatest(Result r1, Result r2)
{
if (r1.dis > r2.dis) return r1;
return r2;
}
}

public class Point3dIndex : Point3d
{
public Point3dIndex(double x, double y, double z, int Idx)
: base(x, y, z)
{
Index = Idx;
}

public int Index;

}

public class ResearchMin
{
List<Point3d> allpoints;
int SplitFactor;
double LatoCubo;

public ResearchMin(List<Point3d> input, double LCubo)
{
allpoints = input;
LatoCubo = LCubo;
}

List<Point3dIndex>[, ,] Zoner;

public Result Research()
{
int ln = allpoints.Count;
if (ln == 1)
{
return new Result(0, 0, 0);
}
if (ln == 2)
{
double dis = allpoints[1] - allpoints[2];
Result rr = new Result(dis, 0, 1);
}

Result res = new Result(double.MaxValue);
double SplitFactorD = (int)(2 * Math.Sqrt(Math.Sqrt(allpoints.Count)));
int SplitFactor = (int)SplitFactorD;

RsMin(SplitFactor, res, false);
RsMin(SplitFactor, res, true);

if (res.dis > (LatoCubo / SplitFactorD))
throw new Exception("Minicubi troppo piccoli (Pochi punti, soluzione classica)");

return res;
}

private Result RsMin(int sp, Result res, bool Sfasa)
{
SplitFactor = sp;

double DimCella = LatoCubo / (double)SplitFactor;

// Initialize Cluster
int SpFSf = SplitFactor + (Sfasa ? 1 : 0);
Zoner = new List<Point3dIndex>[SpFSf, SpFSf, SpFSf];
for (int x = 0; x < SpFSf; x++)
{
for (int y = 0; y < SpFSf; y++)
{
for (int z = 0; z < SpFSf; z++)
{
Zoner[x, y, z] = new List<Point3dIndex>();
}
}
}

double sfasa = (Sfasa) ? (DimCella / 2) : 0;
// Clusterize
int ln = allpoints.Count;
for (int t = 0; t < ln; t++)
{
Point3d p3d = allpoints[t];
Point3dIndex p3di = new Point3dIndex(p3d.x, p3d.y, p3d.z, t);

int nx = Normalize(p3d.x + sfasa, LatoCubo, SplitFactor);
int ny = Normalize(p3d.y + sfasa, LatoCubo, SplitFactor);
int nz = Normalize(p3d.z + sfasa, LatoCubo, SplitFactor);

Zoner[nx, ny, nz].Add(p3di);
}

for (int x = 0; x < SpFSf; x++)
{
for (int y = 0; y < SpFSf; y++)
{
for (int z = 0; z < SpFSf; z++)
{
MinWithinList1(Zoner[x, y, z], res);
}
}
}

return res;
}

private void MinWithinList1(List<Point3dIndex> src, Result CurrentResult)
{
int ln = src.Count;
for (int t = 0; t < ln - 1; t++)
{
Point3dIndex p3di1 = src[t];
for (int u = t + 1; u < ln; u++)
{
Point3dIndex p3di2 = src[u];
double dist = p3di2 - p3di1;

if (dist < CurrentResult.dis)
{
CurrentResult.dis = dist;
CurrentResult.index1 = p3di1.Index;
CurrentResult.index2 = p3di2.Index;
}
}
}
}

private int Normalize(double val, double max, int SplitFactor)
{
return (int)((val / max) * (double)SplitFactor);
}
}

public class ResearchMax2
{
private class rm2
{
public rm2(Point3d rif)
{
riferimento = rif;
}
public Point3dIndex punto;
public Point3d riferimento;
public double cdist = double.MaxValue;
}

List<Point3d> allpoints;
double LatoCubo;

public ResearchMax2(List<Point3d> ap, double lc)
{
allpoints = ap;
LatoCubo = lc;
}

rm2 cmmm = new rm2(new Point3d(0, 0, 0));
rm2 cmmp = new rm2(new Point3d(0, 0, 1));
rm2 cmpm = new rm2(new Point3d(0, 1, 0));
rm2 cmpp = new rm2(new Point3d(0, 1, 1));
rm2 cpmm = new rm2(new Point3d(1, 0, 0));
rm2 cpmp = new rm2(new Point3d(1, 0, 1));
rm2 cppm = new rm2(new Point3d(1, 1, 0));
rm2 cppp = new rm2(new Point3d(1, 1, 1));

public Result Research()
{
double LMez = LatoCubo / 2;
int ln = allpoints.Count;
for (int i = 0; i < ln; i++)
{
Point3d p3d = allpoints[i];
double x = p3d.x;
double y = p3d.y;
double z = p3d.z;
rm2 pref = null;
if (x < LMez && y < LMez && z < LMez) pref = cmmm;
else if (x < LMez && y < LMez && z >= LMez) pref = cmmp;
else if (x < LMez && y >= LMez && z < LMez) pref = cmpm;
else if (x < LMez && y >= LMez && z >= LMez) pref = cmpp;
else if (x >= LMez && y < LMez && z < LMez) pref = cpmm;
else if (x >= LMez && y < LMez && z >= LMez) pref = cpmp;
else if (x >= LMez && y >= LMez && z < LMez) pref = cppm;
else if (x >= LMez && y >= LMez && z >= LMez) pref = cppp;

double dst = p3d - pref.riferimento;
if (dst < pref.cdist)
{
pref.cdist = dst;
pref.punto = new Point3dIndex(p3d.x, p3d.y, p3d.z, i);
}
}

double d1 = cmmm.punto - cppp.punto;
Result r1 = new Result(d1, cmmm.punto.Index, cppp.punto.Index);

double d2 = cmmp.punto - cppm.punto;
Result r2 = new Result(d2, cmmp.punto.Index, cppm.punto.Index);

double d3 = cmpm.punto - cpmp.punto;
Result r3 = new Result(d3, cmpm.punto.Index, cpmp.punto.Index);

double d4 = cmpp.punto - cpmm.punto;
Result r4 = new Result(d4, cmpp.punto.Index, cpmm.punto.Index);

Result CurrentResult = Result.Greatest(r1, r2);
CurrentResult = Result.Greatest(CurrentResult, r3);
CurrentResult = Result.Greatest(CurrentResult, r4);

return CurrentResult;
}
}

public class Universe
{
public List<Point3d> Allpoints;

public Universe()
{
Allpoints = new List<Point3d>();
}

public void Load(string fname)
{
StreamReader sr = new StreamReader(fname);
string stt = sr.ReadLine();
string[] first = stt.Split(new char[] { ' ' }, StringSplitOptions.RemoveEmptyEntries);
int howmany = int.Parse(first[0]);
LatoCubo = double.Parse(first[1]);
string line;
while ((line = sr.ReadLine()) != null)
{
string[] coords = line.Split(new char[] { ' ' }, StringSplitOptions.RemoveEmptyEntries);
Allpoints.Add(new Point3d(double.Parse(coords[0]), double.Parse(coords[1]), double.Parse(coords[2])));
}
sr.Close();
}

public void Save(string fname)
{
StreamWriter sw = new StreamWriter(fname);
sw.WriteLine("{0} {1}", Allpoints.Count, LatoCubo);
foreach (Point3d p3d in Allpoints)
{
sw.WriteLine(p3d.ToString());
}
sw.Close();
}

public double LatoCubo;
public void Generate(int npoints, double dim)
{
LatoCubo = dim;
Allpoints.Clear();
for (int t = 0; t < npoints; t++)
{
Allpoints.Add(Point3d.Rand(dim));
}
}
}
}

[Vincenzo1968]

Rimetto qui il codice, ma dovrebbe essere uguale a prima
...


Niente da fare. Ho sempre lo stesso problema :cry:

:bimbo:

[gugoXX]

Niente da fare. Ho sempre lo stesso problema :cry:

:bimbo:

giuro che a me funziona e non ho mai avuto neppure un problema.
Magari e' un probelma di memoria (non capisco come pero' a me possa funzionare sia con 100.000 che con 1.000.000 di punti).
Posso provare a dimezzare la memoria occupata, ma non so se il problema sta proprio li'

[Vincenzo1968]

giuro che a me funziona e non ho mai avuto neppure un problema.
Magari e' un probelma di memoria (non capisco come pero' a me possa funzionare sia con 100.000 che con 1.000.000 di punti).
Posso provare a dimezzare la memoria occupata, ma non so se il problema sta proprio li'

Se è un problema di memoria, può essere dovuto al fatto che dispongo meno di 1 GB?
Questa è la mia macchina:

AMD Athlon(tm) 64 X2
Dual Core Processor 4800+
2.50 GHz
896 MB di RAM

Microsoft Windows XP Professional (32 bit)
Service Pack 3


Ti posto il progetto Visual Studio e, a parte, l'eseguibile:

Progetto Visual Studio 2008:
http://www.guidealgoritmi.it/public/contest08AGugo.zip

Eseguibile:
http://www.guidealgoritmi.it/public/contest08AGugoBin.zip

[gugoXX]

Se è un problema di memoria, può essere dovuto al fatto che dispongo meno di 1 GB?


Un regalino di Natale? :p
Scrivo io la letterina?

Progetto Visual Studio 2008:
http://www.guidealgoritmi.it/public/contest08AGugo.zip

Eseguibile:
http://www.guidealgoritmi.it/public/contest08AGugoBin.zip

Lo zip dell'eseguibile e' corrotto (problemi di memoria? :D )

[banryu79]

Lo zip dell'eseguibile e' corrotto (problemi di memoria? :D )

No, semplicemente è italiano.

(battuta orribile e qualunquista, lo so, perdonatemi ma oggi in ufficio c'è la desolazione più totale e io sono in modalità cazzeggio-prenatalizio).

[Vincenzo1968]

Un regalino di Natale? :p
Scrivo io la letterina?



Lo zip dell'eseguibile e' corrotto (problemi di memoria? :D )

Li mortan guerieri!!!
L'ho ricompilato più de na vorta, ma mme dà sempre li stessi problemi. Li mortan de li mortan guerieri!!!
Me poi posta l'eseguibile tuo?

:bimbo:

[Vincenzo1968]

Ho ricompilato. Prova con questo:

http://www.guidealgoritmi.it/public/contest08aGugo.zip

:bimbo:

[Tommo]

M'è venuta un'idea per trovare la distanza massima, che dovrebbe essere molto più veloce... non so però se sia esatta:

Invece che provare tutti i "miniscatoli" (che sarebbero octree cells ma vabè :asd:) è facile capire quale è quello che contiene la coppia a distanza massima.

Preso un punto P qualsiasi:
-dividi la "scatola universo" in 8 "sottovolumi" per P
-sono definiti dai 3 piani allineati agli assi passanti per P
-la celletta che contiene il punto più distante è quella contenuta nel sottovolume di volume maggiore, e che a sua volta contiene uno dei vertici del cubo universo.
-Si fa quindi la distanza tra P ed ogni punto contenuto nella celletta trovata, in media meno di 10 iterazioni.

In questa maniera il tutto sarebbe praticamente lineare, e dovrebbe metterci infinitamente meno di ora... però non so come comportarmi nel caso che la cella trovata sia vuota. :stordita:

Cmq per trovare la coppia a distanza massima basta calcolarsi i punti nelle 4 coppie di celle ai lati opposti della diagonale del cubo, credo.

[gugoXX]

M'è venuta un'idea per trovare la distanza massima, che dovrebbe essere molto più veloce... non so però se sia esatta:

Invece che provare tutti i "miniscatoli" (che sarebbero octree cells ma vabè :asd:) è facile capire quale è quello che contiene la coppia a distanza massima.

Preso un punto P qualsiasi:
-dividi la "scatola universo" in 8 "sottovolumi" per P
-sono definiti dai 3 piani allineati agli assi passanti per P
-la celletta che contiene il punto più distante è quella contenuta nel sottovolume di volume maggiore, e che a sua volta contiene uno dei vertici del cubo universo.
-Si fa quindi la distanza tra P ed ogni punto contenuto nella celletta trovata, in media meno di 10 iterazioni.

In questa maniera il tutto sarebbe praticamente lineare, e dovrebbe metterci infinitamente meno di ora... però non so come comportarmi nel caso che la cella trovata sia vuota. :stordita:

Cmq per trovare la coppia a distanza massima basta calcolarsi i punti nelle 4 coppie di celle ai lati opposti della diagonale del cubo, credo.

E' l'algoritmo implementato da me...
Ma e' greedy, nel senso che puo' sbagliare. Il piu' delle volte trova una soluzione molto buona, su grandi numeri (casuali) trova spesso la soluzione ottima.
Ma appunto non c'e' garanzia, puo' anche sbagliare.

[gugoXX]

Ho ricompilato. Prova con questo:

http://www.guidealgoritmi.it/public/contest08aGugo.zip

:bimbo:

Funziona subito...
:nono:

[Vincenzo1968]

Funziona subito...
:nono:

E allora è un problema di memoria? :confused:

[marco.r]

La soluzione di Marco.r, invece, va in crash.

:bimbo:

:mbe: uhm, daro' una occhiata, mi sembra strano perche' non faccio gestione esplicita della memoria. Riesci eventualmente a dirmi se ti va in crash con la ricerca del massimo o del minimo (immagino minimo che e' piu' lenta e usa piu' ram). Puo' essere semplicemente che usa troppa memoria... che dataset hai usato ?

[Vincenzo1968]

:mbe: uhm, daro' una occhiata, mi sembra strano perche' non faccio gestione esplicita della memoria. Riesci eventualmente a dirmi se ti va in crash con la ricerca del massimo o del minimo (immagino minimo che e' piu' lenta e usa piu' ram). Puo' essere semplicemente che usa troppa memoria... che dataset hai usato ?

Ehm,

non saprei. Il programma lavora per un po' e poi va in crash(con il classico messaggio che chiede alll'utente se si vuole inviare l'errore a Microsoft).
Faccio qualche altra prova e ti faccio sapere.

Che dataset ho usato? boh :confused: (se ti riferisci al file ho usato quello di Gugo da 100.000 records).

[Vincenzo1968]

:mbe: uhm, daro' una occhiata, mi sembra strano perche' non faccio gestione esplicita della memoria. Riesci eventualmente a dirmi se ti va in crash con la ricerca del massimo o del minimo (immagino minimo che e' piu' lenta e usa piu' ram). Puo' essere semplicemente che usa troppa memoria... che dataset hai usato ?

Commentando, nel main, la chiamata alla funzione solve_max, non ho nessun problema:

File piccolo:
http://www.guidealgoritmi.it/images/ImgForums/Contest08ASoluzMarcor.jpg

File grande:
http://www.guidealgoritmi.it/images/ImgForums/Contest08ASoluzMarcor2.jpg

[marco.r]

M'è venuta un'idea per trovare la distanza massima, che dovrebbe essere molto più veloce... non so però se sia esatta:

Invece che provare tutti i "miniscatoli" (che sarebbero octree cells ma vabè :asd:) è facile capire quale è quello che contiene la coppia a distanza massima.

Preso un punto P qualsiasi:
-dividi la "scatola universo" in 8 "sottovolumi" per P
-sono definiti dai 3 piani allineati agli assi passanti per P
-la celletta che contiene il punto più distante è quella contenuta nel sottovolume di volume maggiore, e che a sua volta contiene uno dei vertici del cubo universo.
-Si fa quindi la distanza tra P ed ogni punto contenuto nella celletta trovata, in media meno di 10 iterazioni.

In questa maniera il tutto sarebbe praticamente lineare, e dovrebbe metterci infinitamente meno di ora... però non so come comportarmi nel caso che la cella trovata sia vuota. :stordita:

Cmq per trovare la coppia a distanza massima basta calcolarsi i punti nelle 4 coppie di celle ai lati opposti della diagonale del cubo, credo.

Se vuoi rendere il tutto "praticamente lineare", l'unica soluzione che mi viene in mente e' quella di stabilire una euristica che ti permetta di escludere molto velocemente (in tempo costante) la maggior parte dei punti, aggiungendo pochi punti al set dei "papabili" e togliendone il piu' possibile man mano che procedi.
La mia soluzione e' un esempio (se la distanza max trovata e' d, posso escludere un punto che si trovi a meno di 'd' dal piu' lontano angolo in quanto non posso certo trovare punti al di fuori della sfera di quelal dimensione)ma probabilmente se ne possono pensare di altre.

[Vincenzo1968]

Ahò, li mortan guerieri!!!

So' riuscito a superà i tempi de repne ppure pe la ricerca daa distanza massima:

File piccolo:
http://www.guidealgoritmi.it/images/ImgForums/Contest08ASoluzFinale1.jpg

File Grosso:
http://www.guidealgoritmi.it/images/ImgForums/Contest08ASoluzFinale2.jpg

La mia macchina:

AMD Athlon(tm) 64 X2
Dual Core Processor 4800+
2.50 GHz
896 MB di RAM

Microsoft Windows XP Professional (32 bit)
Service Pack 3



El còdigo fuente:

#include <stdio.h>
#include <math.h>
#include <time.h>
#include <malloc.h>

#define BUFFER_SIZE 4096
#define MAX_STACK 100

typedef double ElementType;
typedef ElementType ItemType[3];

typedef struct tagPunti
{
ItemType data;
int index;
} Punti;

typedef struct tagRes
{
Punti P1;
Punti P2;
double Distanza;
int index1;
int index2;
} Res;

typedef struct tagArray
{
int dimensione;
double dmax;
Punti *m;
} Array;

typedef struct tagTree
{
Punti P;
//struct tagTree *father;
struct tagTree *left;
struct tagTree *right;
} Tree;

Tree *NewNode(ItemType key, int index);
Tree *InsertNode(Tree *node, ItemType key, int index);
void FreeTree(Tree *head);


typedef enum tagStati
{
S_ERROR = -1, S0 = 0, S1, S2, S3, S4, S5, S6, S7, S8, S9
} Stati;


// La macro seguente restituisce il numero degli elementi di un array.
// Per esempio, data la seguente dichiarazione:
//
// int array[] = {8, 5, 13, 2, 1};
//
// ARRAY_SIZE(array) restituisce 5

#define ARRAY_SIZE(Array) (sizeof(Array) / sizeof((Array)[0]))

Stati DFA(const char *szFileName, Array *pArray)
{
Stati stato = S0;
FILE *fp;
unsigned char buffer[BUFFER_SIZE];
int numblocks;
int numread;
unsigned char c;
int k, x, j;
int riga, colonna;
char szNum[256];
double num;

unsigned char byteCR = 0xD; // Carriage Return
unsigned char byteLF = 0xA; // Line Feed

fp = fopen(szFileName, "rb");
if ( fp == NULL )
{
printf("Errore nell'apertura del file %s\n", szFileName);
return S_ERROR;
}

if ( fseek(fp, 0, SEEK_END) )
return 0;

numblocks = ftell(fp)/BUFFER_SIZE;
if ( numblocks == 0 )
{
numblocks = 1;
}
else
{
if ( ftell(fp) % BUFFER_SIZE != 0 )
numblocks++;
}

fseek(fp, 0, SEEK_SET);
numread = fread(buffer, 1, BUFFER_SIZE, fp);
if (numread == 0)
{
fclose(fp);
printf("Errore 1 nella lettura del file %s\n", szFileName);
return S_ERROR;
}

pArray->dimensione = 0;
k = 0;
c = *(buffer + k++);
while ( 1 )
{
if ( c >= '0' && c <= '9' )
{
pArray->dimensione = pArray->dimensione * 10 + c - '0';
}
else if ( c == ' ' || c == '\t' )
{
j = 0;
while ( c != byteLF )
{
c = *(buffer + k++);
if ( c >= '0' && c <= '9' )
szNum[j++] = c;
}
szNum[j] = '\0';
pArray->dmax = atof(szNum);
break;
}
else if ( c == '\n' )
{
break;
}

c = *(buffer + k++);
}

pArray->m = (Punti*)malloc(pArray->dimensione*sizeof(Punti));
if ( !pArray->m )
{
printf("Memoria insufficiente.\n");
fclose(fp);
return S_ERROR;
}

riga = colonna = 0;
num = 0;
x = j = 0;
while ( x < numblocks )
{
c = *(buffer + k++);

if ( c == byteLF )
{
pArray->m[riga].index = riga;
riga++;
colonna = 0;
}

switch (stato)
{
case S0:
j = 0;
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S2;
}
else if ( c == '.' )
{
szNum[j++] = c;
stato = S3;
}
else if ( c == '-' || c == '+' )
{
szNum[j++] = c;
stato = S1;
}
else if ( c == ' ' || c == '\t' )
{
while ( c == ' ' || c == '\t' )
{
c = *(buffer + k++);
if ( k >= numread )
break;
}
k--;
}
else if ( c == byteCR || c == byteLF )
{
// nada
}
else
{
printf("Errore: stato S0; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S1:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S2;
}
else if ( c == '.' )
{
szNum[j++] = c;
stato = S3;
}
else
{
printf("Errore: stato S1; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S2:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
}
else if ( c == '.' )
{
szNum[j++] = c;
stato = S4;
}
else if ( c == 'E' || c == 'e' )
{
szNum[j++] = c;
stato = S5;
}
else if ( c == ' ' || c == '\t' || c == byteCR || c == byteLF )
{
szNum[j] = '\0';
pArray->m[riga].data[colonna++] = atof(szNum);
stato = S0;
}
else
{
printf("Errore: stato S2; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S3:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S6;
}
else
{
printf("Errore: stato S3; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S4:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S6;
}
else if ( c == 'E' || c == 'e' )
{
szNum[j++] = c;
stato = S7;
}
else if ( c == ' ' || c == '\t' || c == byteCR || c == byteLF )
{
szNum[j] = '\0';
pArray->m[riga].data[colonna++] = atof(szNum);
stato = S0;
}
else
{
printf("Errore: stato S4; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S5:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S9;
}
else if ( c == '-' || c == '+' )
{
szNum[j++] = c;
stato = S8;
}
else
{
printf("Errore: stato S5; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S6:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S6;
}
else if ( c == 'E' || c == 'e' )
{
szNum[j++] = c;
stato = S7;
}
else if ( c == ' ' || c == '\t' || c == byteCR || c == byteLF )
{
szNum[j] = '\0';
pArray->m[riga].data[colonna++] = atof(szNum);
stato = S0;
}
else
{
printf("Errore: stato S6; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S7:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S9;
}
else if ( c == '-' || c == '+' )
{
szNum[j++] = c;
stato = S8;
}
else
{
printf("Errore: stato S7; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S8:
case S9:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S9;
}
else if ( c == ' ' || c == '\t' || c == byteCR || c == byteLF )
{
szNum[j] = '\0';
pArray->m[riga].data[colonna++] = atof(szNum);
stato = S0;
}
else
{
printf("Errore: stato S9; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;
}

if ( k >= numread )
{
numread = fread(buffer, 1, BUFFER_SIZE, fp);
k = 0;
x++;
}
}

fclose(fp);

return stato;
}

Tree *NewNode(ItemType key, int index)
{
Tree *r;

r = (Tree *) malloc(sizeof(Tree));

if(!r)
{
printf("Memoria insufficiente.\n");
return NULL;
}

r->P.data[0] = key[0];
r->P.data[1] = key[1];
r->P.data[2] = key[2];
r->P.index = index;
//r->father = NULL;
r->left = NULL;
r->right = NULL;

return r;
}

Tree *InsertNode(Tree *node, ItemType key, int index)
{
Tree *pRadice = NULL;
int Level = 0;

if( !node )
{
node = NewNode(key, index);
return node;
}

pRadice = node;

while( 1 )
{
if ( key[Level] < node->P.data[Level] )
{
if ( !node->left )
{
node->left = NewNode(key, index);
//node->left->father = node;
break;
}
node = node->left;
}
else
{
if ( !node->right )
{
node->right = NewNode(key, index);
//node->right->father = node;
break;
}
node = node->right;
}
Level++;
if ( Level > 2 )
Level = 0;
}

node = pRadice;

return node;
}

void FreeTree(Tree *head)
{
Tree *temp1, *temp2;

Tree *stack[MAX_STACK];
int top;

top = 0;

if ( !head )
return;

temp1 = temp2 = head;

while ( temp1 != NULL )
{
for(; temp1->left != NULL; temp1 = temp1->left)
stack[top++] = temp1;

while ( (temp1 != NULL) && (temp1->right == NULL || temp1->right == temp2) )
{
temp2 = temp1;
free(temp2);
if ( top == 0 )
return;
temp1 = stack[--top];
}

stack[top++] = temp1;
temp1 = temp1->right;
}
}

void TreeToArrayMin(Tree *head, Array *pArray)
{
Tree *temp;

Tree *stack[MAX_STACK];
int top = 0;
int k = 0;

if ( !head )
return;

temp = head;

stack[top++] = temp; // Push

while ( top != 0 )
{
temp = stack[--top]; // Pop

//printf("%d -> (%lf,%lf,%lf)\n", temp->index, temp->data[0], temp->data[1], temp->data[2]);
pArray->m[k++] = temp->P;

if ( temp->right != NULL )
stack[top++] = temp->right;
if ( temp->left != NULL )
stack[top++] = temp->left;
}
}

void CalcolaMin(Punti P[], int dim, Res *pMin)
{
int x;
double distanza, distanzaMin;
double diffA, diffB, diffC;

distanzaMin = -1;

for ( x = 1; x < dim; x++ )
{
diffA = P[x-1].data[0] - P[x].data[0];
diffB = P[x-1].data[1] - P[x].data[1];
diffC = P[x-1].data[2] - P[x].data[2];
distanza = diffA*diffA + diffB*diffB + diffC*diffC;
if ( distanza < distanzaMin || distanzaMin == -1 )
{
distanzaMin = distanza;
pMin->P1 = P[x-1];
pMin->P2 = P[x];
pMin->index1 = P[x-1].index;
pMin->index2 = P[x].index;
}
}

pMin->Distanza = sqrt( distanzaMin );
}

void FindMaximumDistance(Array *pArray, Res *pMax)
{
int k, y, x, j;
ElementType min, maxd;
Punti a[8];
Punti r[8];
ElementType diffA, diffB, diffC;
ElementType distanza;

a[0].data[0] = 0;
a[0].data[1] = 0;
a[0].data[2] = 0;

a[1].data[0] = 0;
a[1].data[1] = pArray->dmax;
a[1].data[2] = 0;

a[2].data[0] = 0;
a[2].data[1] = 0;
a[2].data[2] = pArray->dmax;

a[3].data[0] = pArray->dmax;
a[3].data[1] = 0;
a[3].data[2] = 0;

a[4].data[0] = pArray->dmax;
a[4].data[1] = pArray->dmax;
a[4].data[2] = pArray->dmax;

a[5].data[0] = pArray->dmax;
a[5].data[1] = pArray->dmax;
a[5].data[2] = 0;

a[6].data[0] = pArray->dmax;
a[6].data[1] = 0;
a[6].data[2] = pArray->dmax;

a[7].data[0] = 0;
a[7].data[1] = pArray->dmax;
a[7].data[2] = pArray->dmax;

for (y = 0; y < 8; y++)
{
j = x = 0;
min = 1000000000;
for ( k = 0; k < pArray->dimensione; k++ )
{
diffA = pArray->m[k].data[0] - a[y].data[0];
diffB = pArray->m[k].data[1] - a[y].data[1];
diffC = pArray->m[k].data[2] - a[y].data[2];
distanza = diffA*diffA + diffB*diffB + diffC*diffC;

if ( distanza < min )
{
r[y] = pArray->m[k];
min = distanza;
}
}
}

x = y = 0;
maxd = -1000000000;
while( x < 7 )
{
y = x + 1;
while( y < 8 )
{
//distanza = sqrt( pow(P[x].A - P[y].A, 2.0) + pow(P[x].B - P[y].B, 2.0) + pow(P[x].C - P[y].C, 2.0) );
diffA = r[x].data[0] - r[y].data[0];
diffB = r[x].data[1] - r[y].data[1];
diffC = r[x].data[2] - r[y].data[2];
distanza = diffA*diffA + diffB*diffB + diffC*diffC;

if ( distanza > maxd )
{
maxd = distanza;

pMax->P1 = r[x];
pMax->P2 = r[y];

pMax->index1 = r[x].index;
pMax->index2 = r[y].index;
}
++y;
}
++x;
}

pMax->Distanza = sqrt( maxd );
}

int main()
{
Stati stato;
Array myArray;
Res r1, r2;
Tree *pTree = NULL;
int k;
double t1, t2;

clock_t c_start, c_end;

//char *szFileName = "C:\\Temp\\Coords\\Coords.dat";
char *szFileName = "C:\\Temp\\Coords2\\Coords.dat";

c_start = clock();
stato = DFA(szFileName, &myArray);
c_end = clock();
printf("\nTempo impiegato per caricamento dati -> %5.5f secondi\n", (double)(c_end - c_start) / CLOCKS_PER_SEC);

c_start = clock();
for(k = 0; k < myArray.dimensione; k++)
pTree = InsertNode(pTree, myArray.m[k].data, myArray.m[k].index);
TreeToArrayMin(pTree, &myArray);
CalcolaMin(myArray.m, myArray.dimensione, &r1);
c_end = clock();
t1 = (double)(c_end - c_start) / CLOCKS_PER_SEC;

c_start = clock();
FindMaximumDistance(&myArray, &r2);
c_end = clock();
t2 = (double)(c_end - c_start) / CLOCKS_PER_SEC;

printf("\nDistanza Min: P[%d] - P[%d] -> %.15lf\n", r1.index1, r1.index2, r1.Distanza);
printf("\nDistanza Max: P[%d] - P[%d] -> %.15lf\n", r2.index1, r2.index2, r2.Distanza);
printf("\nTempo impiegato per la ricerca della distanza minima -> %5.5f secondi\n", t1);
printf("\nTempo impiegato per la ricerca della distanza massima -> %5.5f secondi\n", t2);

free(myArray.m);
FreeTree(pTree);

return 0;
}


:bimbo:

[Vincenzo1968]

Ho ottimizzato il codice per la ricerca della distanza massima:

Tempi sul file grosso:
http://www.guidealgoritmi.it/images/ImgForums/Contest08ASoluzFinale3.jpg


#include <stdio.h>
#include <math.h>
#include <time.h>
#include <malloc.h>

#define BUFFER_SIZE 4096
#define MAX_STACK 100

typedef double ElementType;
typedef ElementType ItemType[3];

typedef struct tagPunti
{
ItemType data;
int index;
} Punti;

typedef struct tagRes
{
Punti P1;
Punti P2;
double Distanza;
int index1;
int index2;
} Res;

typedef struct tagArray
{
int dimensione;
double dmax;
Punti *m;
} Array;

typedef struct tagTree
{
Punti P;
//struct tagTree *father;
struct tagTree *left;
struct tagTree *right;
} Tree;

Tree *NewNode(ItemType key, int index);
Tree *InsertNode(Tree *node, ItemType key, int index);
void FreeTree(Tree *head);


typedef enum tagStati
{
S_ERROR = -1, S0 = 0, S1, S2, S3, S4, S5, S6, S7, S8, S9
} Stati;


// La macro seguente restituisce il numero degli elementi di un array.
// Per esempio, data la seguente dichiarazione:
//
// int array[] = {8, 5, 13, 2, 1};
//
// ARRAY_SIZE(array) restituisce 5

#define ARRAY_SIZE(Array) (sizeof(Array) / sizeof((Array)[0]))

Stati DFA(const char *szFileName, Array *pArray)
{
Stati stato = S0;
FILE *fp;
unsigned char buffer[BUFFER_SIZE];
int numblocks;
int numread;
unsigned char c;
int k, x, j;
int riga, colonna;
char szNum[256];
double num;

unsigned char byteCR = 0xD; // Carriage Return
unsigned char byteLF = 0xA; // Line Feed

fp = fopen(szFileName, "rb");
if ( fp == NULL )
{
printf("Errore nell'apertura del file %s\n", szFileName);
return S_ERROR;
}

if ( fseek(fp, 0, SEEK_END) )
return 0;

numblocks = ftell(fp)/BUFFER_SIZE;
if ( numblocks == 0 )
{
numblocks = 1;
}
else
{
if ( ftell(fp) % BUFFER_SIZE != 0 )
numblocks++;
}

fseek(fp, 0, SEEK_SET);
numread = fread(buffer, 1, BUFFER_SIZE, fp);
if (numread == 0)
{
fclose(fp);
printf("Errore 1 nella lettura del file %s\n", szFileName);
return S_ERROR;
}

pArray->dimensione = 0;
k = 0;
c = *(buffer + k++);
while ( 1 )
{
if ( c >= '0' && c <= '9' )
{
pArray->dimensione = pArray->dimensione * 10 + c - '0';
}
else if ( c == ' ' || c == '\t' )
{
j = 0;
while ( c != byteLF )
{
c = *(buffer + k++);
if ( c >= '0' && c <= '9' )
szNum[j++] = c;
}
szNum[j] = '\0';
pArray->dmax = atof(szNum);
break;
}
else if ( c == '\n' )
{
break;
}

c = *(buffer + k++);
}

pArray->m = (Punti*)malloc(pArray->dimensione*sizeof(Punti));
if ( !pArray->m )
{
printf("Memoria insufficiente.\n");
fclose(fp);
return S_ERROR;
}

riga = colonna = 0;
num = 0;
x = j = 0;
while ( x < numblocks )
{
c = *(buffer + k++);

if ( c == byteLF )
{
pArray->m[riga].index = riga;
riga++;
colonna = 0;
}

switch (stato)
{
case S0:
j = 0;
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S2;
}
else if ( c == '.' )
{
szNum[j++] = c;
stato = S3;
}
else if ( c == '-' || c == '+' )
{
szNum[j++] = c;
stato = S1;
}
else if ( c == ' ' || c == '\t' )
{
while ( c == ' ' || c == '\t' )
{
c = *(buffer + k++);
if ( k >= numread )
break;
}
k--;
}
else if ( c == byteCR || c == byteLF )
{
// nada
}
else
{
printf("Errore: stato S0; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S1:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S2;
}
else if ( c == '.' )
{
szNum[j++] = c;
stato = S3;
}
else
{
printf("Errore: stato S1; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S2:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
}
else if ( c == '.' )
{
szNum[j++] = c;
stato = S4;
}
else if ( c == 'E' || c == 'e' )
{
szNum[j++] = c;
stato = S5;
}
else if ( c == ' ' || c == '\t' || c == byteCR || c == byteLF )
{
szNum[j] = '\0';
pArray->m[riga].data[colonna++] = atof(szNum);
stato = S0;
}
else
{
printf("Errore: stato S2; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S3:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S6;
}
else
{
printf("Errore: stato S3; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S4:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S6;
}
else if ( c == 'E' || c == 'e' )
{
szNum[j++] = c;
stato = S7;
}
else if ( c == ' ' || c == '\t' || c == byteCR || c == byteLF )
{
szNum[j] = '\0';
pArray->m[riga].data[colonna++] = atof(szNum);
stato = S0;
}
else
{
printf("Errore: stato S4; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S5:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S9;
}
else if ( c == '-' || c == '+' )
{
szNum[j++] = c;
stato = S8;
}
else
{
printf("Errore: stato S5; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S6:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S6;
}
else if ( c == 'E' || c == 'e' )
{
szNum[j++] = c;
stato = S7;
}
else if ( c == ' ' || c == '\t' || c == byteCR || c == byteLF )
{
szNum[j] = '\0';
pArray->m[riga].data[colonna++] = atof(szNum);
stato = S0;
}
else
{
printf("Errore: stato S6; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S7:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S9;
}
else if ( c == '-' || c == '+' )
{
szNum[j++] = c;
stato = S8;
}
else
{
printf("Errore: stato S7; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;

case S8:
case S9:
if ( c >= '0' && c <= '9' )
{
szNum[j++] = c;
stato = S9;
}
else if ( c == ' ' || c == '\t' || c == byteCR || c == byteLF )
{
szNum[j] = '\0';
pArray->m[riga].data[colonna++] = atof(szNum);
stato = S0;
}
else
{
printf("Errore: stato S9; riga %d; colonna %d;\n", riga, colonna);
fclose(fp);
return S_ERROR;
}
break;
}

if ( k >= numread )
{
numread = fread(buffer, 1, BUFFER_SIZE, fp);
k = 0;
x++;
}
}

fclose(fp);

return stato;
}

Tree *NewNode(ItemType key, int index)
{
Tree *r;

r = (Tree *) malloc(sizeof(Tree));

if(!r)
{
printf("Memoria insufficiente.\n");
return NULL;
}

r->P.data[0] = key[0];
r->P.data[1] = key[1];
r->P.data[2] = key[2];
r->P.index = index;
//r->father = NULL;
r->left = NULL;
r->right = NULL;

return r;
}

Tree *InsertNode(Tree *node, ItemType key, int index)
{
Tree *pRadice = NULL;
int Level = 0;

if( !node )
{
node = NewNode(key, index);
return node;
}

pRadice = node;

while( 1 )
{
if ( key[Level] < node->P.data[Level] )
{
if ( !node->left )
{
node->left = NewNode(key, index);
//node->left->father = node;
break;
}
node = node->left;
}
else
{
if ( !node->right )
{
node->right = NewNode(key, index);
//node->right->father = node;
break;
}
node = node->right;
}
Level++;
if ( Level > 2 )
Level = 0;
}

node = pRadice;

return node;
}

void FreeTree(Tree *head)
{
Tree *temp1, *temp2;

Tree *stack[MAX_STACK];
int top;

top = 0;

if ( !head )
return;

temp1 = temp2 = head;

while ( temp1 != NULL )
{
for(; temp1->left != NULL; temp1 = temp1->left)
stack[top++] = temp1;

while ( (temp1 != NULL) && (temp1->right == NULL || temp1->right == temp2) )
{
temp2 = temp1;
free(temp2);
if ( top == 0 )
return;
temp1 = stack[--top];
}

stack[top++] = temp1;
temp1 = temp1->right;
}
}

void TreeToArrayMin(Tree *head, Array *pArray)
{
Tree *temp;

Tree *stack[MAX_STACK];
int top = 0;
int k = 0;

if ( !head )
return;

temp = head;

stack[top++] = temp; // Push

while ( top != 0 )
{
temp = stack[--top]; // Pop

//printf("%d -> (%lf,%lf,%lf)\n", temp->index, temp->data[0], temp->data[1], temp->data[2]);
pArray->m[k++] = temp->P;

if ( temp->right != NULL )
stack[top++] = temp->right;
if ( temp->left != NULL )
stack[top++] = temp->left;
}
}

void CalcolaMin(Punti P[], int dim, Res *pMin)
{
int x;
double distanza, distanzaMin;
double diffA, diffB, diffC;

distanzaMin = -1;

for ( x = 1; x < dim; x++ )
{
diffA = P[x-1].data[0] - P[x].data[0];
diffB = P[x-1].data[1] - P[x].data[1];
diffC = P[x-1].data[2] - P[x].data[2];
distanza = diffA*diffA + diffB*diffB + diffC*diffC;
if ( distanza < distanzaMin || distanzaMin == -1 )
{
distanzaMin = distanza;
pMin->P1 = P[x-1];
pMin->P2 = P[x];
pMin->index1 = P[x-1].index;
pMin->index2 = P[x].index;
}
}

pMin->Distanza = sqrt( distanzaMin );
}

void FindMaximumDistance(Array *pArray, Res *pMax)
{
int k, y, x;
ElementType maxd;
ElementType min[8];
Punti a[8];
Punti r[8];
ElementType diffA, diffB, diffC;
ElementType distanza;

a[0].data[0] = 0;
a[0].data[1] = 0;
a[0].data[2] = 0;

a[1].data[0] = 0;
a[1].data[1] = pArray->dmax;
a[1].data[2] = 0;

a[2].data[0] = 0;
a[2].data[1] = 0;
a[2].data[2] = pArray->dmax;

a[3].data[0] = pArray->dmax;
a[3].data[1] = 0;
a[3].data[2] = 0;

a[4].data[0] = pArray->dmax;
a[4].data[1] = pArray->dmax;
a[4].data[2] = pArray->dmax;

a[5].data[0] = pArray->dmax;
a[5].data[1] = pArray->dmax;
a[5].data[2] = 0;

a[6].data[0] = pArray->dmax;
a[6].data[1] = 0;
a[6].data[2] = pArray->dmax;

a[7].data[0] = 0;
a[7].data[1] = pArray->dmax;
a[7].data[2] = pArray->dmax;

min[0] = 1000000000;
min[1] = 1000000000;
min[2] = 1000000000;
min[3] = 1000000000;
min[4] = 1000000000;
min[5] = 1000000000;
min[6] = 1000000000;
min[7] = 1000000000;
for ( k = 0; k < pArray->dimensione; k++ )
{
diffA = pArray->m[k].data[0] - a[0].data[0];
diffB = pArray->m[k].data[1] - a[0].data[1];
diffC = pArray->m[k].data[2] - a[0].data[2];
distanza = diffA*diffA + diffB*diffB + diffC*diffC;
if ( distanza < min[0] )
{
r[0] = pArray->m[k];
min[0] = distanza;
}

diffA = pArray->m[k].data[0] - a[1].data[0];
diffB = pArray->m[k].data[1] - a[1].data[1];
diffC = pArray->m[k].data[2] - a[1].data[2];
distanza = diffA*diffA + diffB*diffB + diffC*diffC;
if ( distanza < min[1] )
{
r[1] = pArray->m[k];
min[1] = distanza;
}

diffA = pArray->m[k].data[0] - a[2].data[0];
diffB = pArray->m[k].data[1] - a[2].data[1];
diffC = pArray->m[k].data[2] - a[2].data[2];
distanza = diffA*diffA + diffB*diffB + diffC*diffC;
if ( distanza < min[2] )
{
r[2] = pArray->m[k];
min[2] = distanza;
}

diffA = pArray->m[k].data[0] - a[3].data[0];
diffB = pArray->m[k].data[1] - a[3].data[1];
diffC = pArray->m[k].data[2] - a[3].data[2];
distanza = diffA*diffA + diffB*diffB + diffC*diffC;
if ( distanza < min[3] )
{
r[3] = pArray->m[k];
min[3] = distanza;
}

diffA = pArray->m[k].data[0] - a[4].data[0];
diffB = pArray->m[k].data[1] - a[4].data[1];
diffC = pArray->m[k].data[2] - a[4].data[2];
distanza = diffA*diffA + diffB*diffB + diffC*diffC;
if ( distanza < min[4] )
{
r[4] = pArray->m[k];
min[4] = distanza;
}

diffA = pArray->m[k].data[0] - a[5].data[0];
diffB = pArray->m[k].data[1] - a[5].data[1];
diffC = pArray->m[k].data[2] - a[5].data[2];
distanza = diffA*diffA + diffB*diffB + diffC*diffC;
if ( distanza < min[5] )
{
r[5] = pArray->m[k];
min[5] = distanza;
}

diffA = pArray->m[k].data[0] - a[6].data[0];
diffB = pArray->m[k].data[1] - a[6].data[1];
diffC = pArray->m[k].data[2] - a[6].data[2];
distanza = diffA*diffA + diffB*diffB + diffC*diffC;
if ( distanza < min[6] )
{
r[6] = pArray->m[k];
min[6] = distanza;
}

diffA = pArray->m[k].data[0] - a[7].data[0];
diffB = pArray->m[k].data[1] - a[7].data[1];
diffC = pArray->m[k].data[2] - a[7].data[2];
distanza = diffA*diffA + diffB*diffB + diffC*diffC;
if ( distanza < min[7] )
{
r[7] = pArray->m[k];
min[7] = distanza;
}
}

x = y = 0;
maxd = -1000000000;
while( x < 7 )
{
y = x + 1;
while( y < 8 )
{
//distanza = sqrt( pow(P[x].A - P[y].A, 2.0) + pow(P[x].B - P[y].B, 2.0) + pow(P[x].C - P[y].C, 2.0) );
diffA = r[x].data[0] - r[y].data[0];
diffB = r[x].data[1] - r[y].data[1];
diffC = r[x].data[2] - r[y].data[2];
distanza = diffA*diffA + diffB*diffB + diffC*diffC;

if ( distanza > maxd )
{
maxd = distanza;

pMax->P1 = r[x];
pMax->P2 = r[y];

pMax->index1 = r[x].index;
pMax->index2 = r[y].index;
}
++y;
}
++x;
}

pMax->Distanza = sqrt( maxd );
}

int main()
{
Stati stato;
Array myArray;
Res r1, r2;
Tree *pTree = NULL;
int k;
double t1, t2;

clock_t c_start, c_end;

//char *szFileName = "C:\\Temp\\Coords\\Coords.dat";
char *szFileName = "C:i\\Temp\\Coords2\\Coords.dat";

c_start = clock();
stato = DFA(szFileName, &myArray);
c_end = clock();
printf("\nTempo impiegato per caricamento dati -> %5.5f secondi\n", (double)(c_end - c_start) / CLOCKS_PER_SEC);

c_start = clock();
for(k = 0; k < myArray.dimensione; k++)
pTree = InsertNode(pTree, myArray.m[k].data, myArray.m[k].index);
TreeToArrayMin(pTree, &myArray);
CalcolaMin(myArray.m, myArray.dimensione, &r1);
c_end = clock();
t1 = (double)(c_end - c_start) / CLOCKS_PER_SEC;

c_start = clock();
FindMaximumDistance(&myArray, &r2);
c_end = clock();
t2 = (double)(c_end - c_start) / CLOCKS_PER_SEC;

printf("\nDistanza Min: P[%d] - P[%d] -> %.15lf\n", r1.index1, r1.index2, r1.Distanza);
printf("\nDistanza Max: P[%d] - P[%d] -> %.15lf\n", r2.index1, r2.index2, r2.Distanza);
printf("\nTempo impiegato per la ricerca della distanza minima -> %5.5f secondi\n", t1);
printf("\nTempo impiegato per la ricerca della distanza massima -> %5.5f secondi\n", t2);

free(myArray.m);
FreeTree(pTree);

return 0;
}

[marco.r]

Commentando, nel main, la chiamata alla funzione solve_max, non ho nessun problema:

Ho cambiato la struttura dati. Fammi sapere se ti da ancora problemi



#include <iostream>
#include <vector>
#include <fstream>
#include <set>
#include <iterator>
#include <cmath>
#include <algorithm>
#include <cassert>

using std::cout;
using std::endl;
using std::make_pair;
using std::set;
using std::pair;
using std::vector;

clock_t start;

void start_clock()
{
start = clock();
}

double stop_clock()
{
clock_t new_start = clock();
float result = ((double)(new_start - start))/CLOCKS_PER_SEC;
start = new_start;
return result;
}

struct point
{
point(){}
point(double _x, double _y, double _z ):x(_x),y(_y),z(_z){}
double x,y,z;
};

std::vector<point> data;
typedef std::vector<int> FeasibleSet;
FeasibleSet feasible_set;
//std::set<int> feasible_set;
point corners[8];

double max_distance = 0.0;
std::pair<int,int> max_pair;
double size;
int dim;



std::istream& operator >> (std::istream& in, point& p )
{
in >> p.x >> p.y >> p.z;
return in;
}

std::ostream& operator << (std::ostream& out, const point& p )
{
out << p.x << ' ' << p.y << ' ' << p.z;
return out;
}


double distance(const point& p1, const point& p2 )
{
double dx,dy,dz;
dx = p1.x - p2.x;
dy = p1.y - p2.y;
dz = p1.z - p2.z;
return dx*dx + dy*dy + dz*dz;
}

pair<int,double> most_distant(const point& p)
{
double dist = 0;
int idx = -1;
for ( FeasibleSet::const_iterator i = feasible_set.begin();
i != feasible_set.end();
++i )
{
double d = distance(p,data[*i]);
if ( d > dist )
{
dist = d;
idx = *i;
}
}
return make_pair(idx,dist);
}

bool feasible( const point& p )
{
for ( int i=0 ; i<8 ; ++i )
{
if ( distance(p,corners[i]) > max_distance )
return true;
}
return false;

}


void init( double size )
{
corners[0] = point( 0, 0, 0 );
corners[1] = point( 0, 0, 1 );
corners[2] = point( 0, 1, 0 );
corners[3] = point( 0, 1, 1 );
corners[4] = point( 1, 0, 0 );
corners[5] = point( 1, 0, 1 );
corners[6] = point( 1, 1, 0 );
corners[7] = point( 1, 1, 1 );

}

bool unfeasible( const int& x )
{
return !feasible(data[x]);
}


void keep_element( FeasibleSet& fs, FeasibleSet::iterator& i )
{
++i;
}

void remove_element( FeasibleSet& fs, FeasibleSet::iterator i )
{
*i = fs.back();
fs.pop_back();
}

void remove_unfeasible()
{
FeasibleSet::iterator i = feasible_set.begin();
while ( i != feasible_set.end() )
{
if ( unfeasible(*i) )
{
remove_element( feasible_set, i );
}
else
{
keep_element( feasible_set, i );
}
}
}

void update_feasible_set( const point& p, int index )
{
if ( feasible(p) )
{
pair<int,double> x = most_distant(p);
const int& idx = x.first;
const double& d = x.second;
if ( d > max_distance )
{
max_distance = d;
max_pair = std::make_pair( index, idx );
}

remove_unfeasible();
feasible_set.push_back(index);
}
}

void solve_max()
{
feasible_set.push_back( 0 );
feasible_set.push_back( 1 );
max_distance = distance( data[0], data[1] );
max_pair = make_pair(0,1);
for ( int i=2 ; i<data.size() ; ++i )
{
update_feasible_set( data[i], i );
}
cout << "Maximum distance is " << sqrt(max_distance) << " between [" << max_pair.first << ',' << max_pair.second << "]" << endl;
}


int nblocks( int points )
{
// voglio avere mediamente 3 punti per blocco

return int( pow( points/3.0, 1.0/3.0 ) )+1;
}

float blocksize( int points, float size )
{
return size/nblocks(points);
}

int block_index(int blocknum , int xidx,int yidx, int zidx)
{
const int& result = xidx*blocknum*blocknum + yidx*blocknum + zidx;
return result;
}
int point_block( float size, float blocksize, int blocknum, const point& p )
{
int xidx,yidx,zidx;
xidx = int( p.x * size/blocksize );
yidx = int( p.y * size/blocksize );
zidx = int( p.z * size/blocksize );
return block_index(blocknum,xidx,yidx,zidx);
}


void throw_in_buckets( vector< set<int> >& buckets, const vector<point>& points,
float size, float blocksize, int blocknum )
{
for ( int i=0 ; i<points.size() ; ++i )
{
int n = point_block( size, blocksize, blocknum, points[i] );
buckets[n].insert(i);
}
}

void find_closest_pair( const vector< set<int> >& buckets, int index1,
float& min_distance, int& first_index, int& second_index )
{
const set<int>& from = buckets[index1];

for ( set<int>::const_iterator i = from.begin();
i != from.end();
++i )
{
for ( set<int>::const_iterator j = from.begin();
j != from.end();
++j )
{
if ( *i == *j )
continue;
float d = distance( data[*i], data[*j] );
if ( d < min_distance )
{
first_index = *i;
second_index = *j;
min_distance = d;
}
}
}
}


void solve_min()
{
// totale confronti : N^2
// se divido in K blocchi con in mesia (N/K) valori:
// K * (N/K)^2 confronti = N^2/K
// prendo K = N/100 = N

int blocknum = nblocks(dim);
float bsize = blocksize(dim,size);
int first_index = 0;
int second_index = 1;
float min_distance = distance( data[0], data[1] );

for ( int a=0; a<3 ; ++a )
{
float bucket_delta = 1.0/(a+1) -1;
vector< set<int> > buckets;
buckets.resize( blocknum*blocknum*blocknum );
throw_in_buckets( buckets, data, size, bsize, blocknum );
for ( int i=0 ; i<blocknum-1 ; ++i )
{
for ( int j=0 ; j<blocknum-1; ++j )
{
for ( int k=0; k<blocknum-1; ++k )
{
for ( int n=0 ; n<8; ++n )
{
int index1 = block_index( blocknum,i,j,k );
find_closest_pair( buckets, index1, min_distance, first_index, second_index );
}
}
}
}
}
cout << "Minimum distance is " << sqrt(min_distance) << " between [" << first_index << ',' << second_index << "]" << endl;
}

int main(int argc, char** argv)
{
const char* filename = argc > 1 ? argv[1] : "Coords.dat";

std::ifstream input(filename);
float load_time,min_time,max_time;

start_clock();
input >> dim >> size;
init(size);

data.resize(dim);
std::copy( std::istream_iterator<point>(input),
std::istream_iterator<point>(),
data.begin() );
load_time = stop_clock();
solve_max();
max_time = stop_clock();
solve_min();
min_time = stop_clock();

cout << "Execution time: " << endl;
cout << "Data load " << load_time << endl;
cout << "Minimum " << min_time << endl;
cout << "Maximum " << max_time << endl;
cout << "====================" << endl;
cout << "Total " << load_time + min_time + max_time << endl;
}

[Vincenzo1968]

Ho cambiato la struttura dati. Fammi sapere se ti da ancora problemi


Adesso non dà nessun problema. Questo è l'output sul file grosso:

http://www.guidealgoritmi.it/images/ImgForums/Contest08ASoluzMarcor3.jpg

Ho dovuto aggiungere questo include


#include <ctime>


perché Visual C++ 2008 dava un errore in compilazione.

Grazie mille ;)

[Vincenzo1968]

Continuo ad avere problemi con la versione C#:

http://www.guidealgoritmi.it/images/ImgForums/Contest08ErroreGugo2.jpg

:confused:

[Vincenzo1968]

Gugo ha capito qual è il problema:

Il file e' stato prodotto su un sistema anglosassone, in cui il punto decimale e' appunto un punto decimale.
Invece lo stai leggendo da un sistema italiano, in cui il punto decimale e' invece una virgola... ergo quello che dovrebbe essere 0.123456789 diventa 1234567

e sta aggiustando il codice.
Nel frattempo, ho preso i tempi sostituendo, nei files, il punto con la virgola:

Tempi sul file grosso:
http://www.guidealgoritmi.it/images/ImgForums/Contest08ASoluzGugo.jpg

;)

[Vincenzo1968]

Classifica per la ricerca della distanza minima:
http://www.guidealgoritmi.it/images/ImgForums/Contest08AClassificaMin.jpg

Classifica per la ricerca della distanza massima:
http://www.guidealgoritmi.it/images/ImgForums/Contest08AClassificaMax.jpg

La mia macchina:

AMD Athlon(tm) 64 X2
Dual Core Processor 4800+
2.50 GHz
896 MB di RAM

Microsoft Windows XP Professional (32 bit)
Service Pack 3