*
* KUMAC pour demontrer la methode locale des cluster
*                                      MP 6/95
opt nstat
opt nbox                                                                        
*                              Histogrammes pour la presentation
his/cre/2d 500 'Distribution x vs. y ' 15 -1. 2. 15 -2. 1.
*                              Lire l'echantillon (voir coulomb.kumac)
ne = 500
vec/cre x1([ne])
vec/cre y1([ne])
vec/read x1 yots1.txt
vec/read y1 tots1.txt
*                              Changer les correllation en laissant V1=V2
sigma x2 = x1 + .5
sigma y2 = y1 - .3
vec/hfill y1%x1 500
vec/hfill y2%x2 500
*                            Extraction des seeds (bins les plus popules)
vec/cre gg(15,15) r
vec/cre null(1) r 0
his/get/con 500 gg
sigma nm=lvmax(gg)
sigma nx=int(mod(nm,15))   ;  sigma ny=int((nm-nx)/15+1)
sigma xs1=-1.+(nx-.5)*.2   ;  sigma ys1=-2.+(ny-.5)*.2
sigma print nx(1),ny(1),xs1(1),ys1(1)
vec/copy null(1:1) gg(8,9)
sigma nm=lvmax(gg)
sigma nx=mod(nm,15)        ;  sigma ny=(nm-nx)/15+1
sigma xs2=-1.+(nx-.5)*.2   ;  sigma ys2=-2.+(ny-.5)*.2
sigma print nx(1),ny(1),xs2(1),ys2(1)
*                              Presentation echantillon et seed
zone 1 1
his/plot 500 'box'
text xs1 ys1 's' .5 ! c    ;  text xs2 ys2 's' .5 ! c
wait
*                            Mettre deux classes dans le meme vecteur
vec/cre xx(1000) r         ;  vec/cre yy(1000) r
vec/copy x1(1:500) xx(1:500)
vec/copy y1(1:500) yy(1:500)
vec/copy x2(1:500) xx(501:1000)
vec/copy y2(1:500) yy(501:1000)
*                              Trouver cluster 1 a partir du seed
*                              V inconnu: coupure en distance brute
sigma dr = sqrt((xx-xs1)**2+(yy-ys1)**2)
*  sigma dr = sqrt(65.*(xx-xs1)**2-2.*64.*(xx-xs1)*(yy-ys1)+85.*(yy-ys1)**2)
sigma ll=lvmin(dr)      ;   sigma dl=dr(ll)
sigma xs1 = xx(ll)      ;   sigma ys1 = yy(ll)
sigma dr=(dr le dr(ll))*1000.
sigma xx=xx+dr
do i=1,600
   sigma dr = sqrt((xx-xs1)**2+(yy-ys1)**2)
*  sigma dr = sqrt(65.*(xx-xs1)**2-2.*64.*(xx-xs1)*(yy-ys1)+85.*(yy-ys1)**2)
   sigma ll=lvmin(dr)   ;   sigma dl=dr(ll)
   if (dl>3.) then
      breakl
   endif
   sigma xs1=(xs1*[i]+xx(ll))/([i]+1)
   sigma ys1=(ys1*[i]+yy(ll))/([i]+1)
   sigma dr=(dr le dr(ll))*1000.
   sigma xx=xx+dr
enddo
*                              Trouver cluster 2 a partir du seed
*                              V connu: coupure  a 3 sigma
sigma dr = sqrt(65.*(xx-xs2)**2-2.*64.*(xx-xs2)*(yy-ys2)+85.*(yy-ys2)**2)
sigma ll=lvmin(dr)      ;   sigma dl=dr(ll)
sigma xs2 = xx(ll)      ;   sigma ys2 = yy(ll)
sigma dr=(dr le dr(ll))*2000.
sigma xx=xx+dr
do i=1,600
   sigma dr = sqrt(65.*(xx-xs2)**2-2.*64.*(xx-xs2)*(yy-ys2)+85.*(yy-ys2)**2)
   sigma ll=lvmin(dr)   ;   sigma dl=dr(ll)
   if (dl>3.) then
      breakl
   endif
   sigma xs2=(xs2*[i]+xx(ll))/([i]+1)
   sigma ys2=(ys2*[i]+yy(ll))/([i]+1)
   sigma dr=(dr le dr(ll))*2000.
   sigma xx=xx+dr
enddo
*                               Presentation cluster
his/plot 500 'box'
set pmci 2   ;   set mtyp 2   ;   pmarker 1000 xx yy
sigma dr= array(1000)*1000.   ;   sigma xx= xx-dr
set pmci 3   ;   set mtyp 3
pmarker 1000 xx yy            ;   text xs1 ys1 '1' .5 ! c
sigma xx= xx-dr
set pmci 4   ;   set mtyp 4
pmarker 1000 xx yy            ;   text xs2 ys2 '2' .5 ! c