[Vari] Contest 8: Compattamento Barionico

[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 (!) () 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]

Originariamente inviato da marco.r

Confronto sempre ogni miniscatolo (!) () 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]

Originariamente inviato da gugoXX

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%

[marco.r]

Originariamente inviato da gugoXX

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

[gugoXX]

Originariamente inviato da 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

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


file da 1.000.000 de punti



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



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):


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

[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;
}

[banryu79]

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

[Vincenzo1968]

Originariamente inviato da banryu79

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

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



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

[^TiGeRShArK^]

Originariamente inviato da Vincenzo1968

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



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

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

[Vincenzo1968]

Originariamente inviato da ^TiGeRShArK^

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

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)