;-First Guess ; ahora 5Sep sobre subsets validos errCutYJ = 0.075 ; 0.25 errCutYH = 0.10 ; 0.35 errCutYK = 0.125 ; 0.35 errCutJH = 0.125 ; 0.25 errCutJK = 0.15 ; 0.25 errCutHK = 0.15 ; 0.15 ;-More stringent ; errCutYJ = 0.1 ; errCutYH = 0.1 ; errCutYK = 0.1 ; errCutJH = 0.1 ; errCutJK = 0.1 ; errCutHK = 0.1 ;-super stupidly noisy *********** ; errCutYJ = 5.0 ; errCutYH = 5.0 ; errCutYK = 5.0 ; errCutJH = 5.0 ; errCutJK = 5.0 ; errCutHK = 5.0 medCut = median( [errCutYJ,errCutYH,errCutYK,errCutJH,errCutJK,errCutHK] ) ;-User choice numSel = 1 ;- cuantos l=colores minimo se tienen que usar: 2, 3 o 4 NONON ; no hay todos los casos referirse al histograma ;---------- Directories spawn, 'hostname', host case strTrim(host) of 'hyperion': dirMOVIS = '/home/bcarry/data/mining/movis/data/' 'endymion': dirMOVIS = '/home/bcarry/data/mining/movis/data/' 'gnanzi': dirMOVIS = '/home/mpajuelo/W/Vista/' 'IMCCE5': dirMOVIS = 'PON_TU_PATH' else: stop endcase ;-What to do showGraphBDM = 0 ;-X createGraphBDM = 0 ;-PNG/EPS showGraphMOVIS = 0 ;-X createGraphMOVIS = 1 ;-PNG/EPS checkTaxo = 0 ;-Classify BDM colors instead of MOVIS (-> auto-check) onlyValidFig = 0 ;-Re run barplot code only - only for checktaxo=1 ; bdmChoice = 'T' ;-Select a single class among BDM (only if checkTaxo eq 1) ;-Plot Parameters ccRang = [ [ 0.10, 0.80], $ ;- Y-J [ 0.35, 1.40], $ ;- Y-H [ 0.45, 1.60], $ ;- Y-K [ 0.00, 0.90], $ ;- J-H [-0.10, 1.10], $ ;- J-K [-0.15, 0.30] ] ;- H-K ccTickInt= [0.2, 0.2, 0.2, 0.2, 0.2, 0.1 ] ccTickMin= [ 4 , 4 , 4 , 4 , 4 , 4 ] if onlyValidFig eq 1 then goto, j2xm ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--- TAG --- I --- Define classes ------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--I.1--Class Labels -------------------------------------------------------------------- ;- 1 2 3 4 5 6 7 8 9 10 11 cLabel=['A','B','C','D','K','L','Q','S','T','X','V'] nbClass = n_elements( cLabel ) print, 'List of classes considered here ('+strtrim(string(nbClass),2)+') ' print, cLabel ;--I.2--Class Look'n'Feel --------------------------------------------------------------- cColor=['Sienna' , $ ;-- A 'Gray' , $ ;-- B 'Charcoal' , $ ;-- C 'Dark Gray' , $ ;-- D 'Green' , $ ;-- K 'Navy', $ ;-- L 'Orange' , $ ;-- Q 'Cornflower blue', $ ;-- S 'PUR8' , $ ;-- T 'Magenta' , $ ;-- X 'Red' , $ ;-- V 'light grey' ] ;-- U ;+++++++++ cSymb =['Filled circle' , $ ;-- A 'Filled circle', $ ;-- B 'Open circle', $ ;-- C 'Filled circle', $ ;-- D 'Open circle', $ ;-- K 'Open square', $ ;-- L 'Filled star', $ ;-- Q 'Open square', $ ;-- S 'Star' , $ ;-- T 'Open diamond', $ ;-- X 'Filled Circle', $ ;-- V 'open circle' ] ;-- U ; ++++++++++ SSymb =[ 0.8 ,$ ;-- A 0.8 ,$ ;-- B 0.8 ,$ ;-- C 0.8 ,$ ;-- D 0.8 ,$ ;-- K 0.8 ,$ ;-- L 1.5 ,$ ;-- Q 0.8 ,$ ;-- S 1.5 ,$ ;-- T 0.8 ,$ ;-- X 0.8 ,$ ;-- V 0.5 ] ;-- U ; ++++++++++ ;--I.X--Useful stuff -------------------------------------------------------------------- theta=2*!PI*findgen(361)/360 ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--- TAG --- II --- Read Datasets ------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--II.1--Read Bus-DeMeo Summary Color --------------------------------------------------- data='(A,A,F,F,F,I,I,F,F,F,I,I,F,F,F,I,I,F,F,F,I,I,F,F,F,I,I,F,F,F,I,I)' readcol, dirMOVIS+'bdm-col.csv',f=data,$ class,ID,YmeJ,dYmeJ,YmeJDt,NsY,NsJ,$ YmeH,dYmeH,YmeHDt,NsY,NsH,$ YmeKs,dYmeKs,YmeKsDt,NsY,NsKs,$ JmeH,dJmeH,JmeHDt,NsJ,NsH,$ JmeKs,dJmeKs,JmeKsDt,NsJ,NsKs,$ HmeKs,dHmeKs,HmeKsDt,NsH,NsKs,$ delimiter=',', /Silent bdmCol=[[YmeJ],[YmeH],[YmeKs],[JmeH],[JmeKs],[HmeKs]] uClass = uniq(class, sort(class)) print, 'List of classes in Bus-DeMeo color file ('+strtrim(string(n_elements(uClass)),2)+') ' print, class[uClass] ;--II.2--Read MOVIS-C csv File ---------------------------------------------------------- if checkTaxo eq 0 then begin ;--II.2.1--Blind read movis = vistaReadMOVIS( dirMOVIS+'MOVIS-C.csv' ) nbIn = n_elements(movis) ;--II.2.2--Selection based on threshold selectMOVIS = where( movis.YmJ.unc le errCutYJ and $ movis.YmH.unc le errCutYH and $ movis.YmK.unc le errCutYK and $ movis.JmH.unc le errCutJH and $ movis.JmK.unc le errCutJK and $ movis.HmK.unc le errCutHK, nbMOVIS) ;--II.2.3--Selection based on threshold movis = movis[selectMOVIS] ssoCol = [[movis.YmJ.val],[movis.YmH.val],[movis.YmK.val],[movis.JmH.val],[movis.JmK.val],[movis.HmK.val]] nbSSO = n_elements( movis ) print, 'Number of asteroids in MOVIS:', nbSSO ssoNum = movis.num ssoName= movis.name ssoKnown= replicate('U',nbSSO) endif else begin ;-read BusDeMeo name readcol, dirMOVIS+'bdm.csv', bdmNum, bdmName, bdmC1, bdmC2, bdmC3, bdmDate, format='(L,A,A,A,A,A)', delimiter=',', /Silent nbIn = n_elements(bdmNum) nbMOVIS = nbIn ssoNum = bdmNum ssoName= bdmName ssoKnown= strtrim(bdmC3,2) ssoCol = bdmCol nbSSO = n_elements( class ) if keyword_set(bdmChoice) then begin select = where( strCmp( bdmChoice, strMid(strTrim(bdmC3,2),0,1) ), nbSSO ) ssoNum = ssoNum[select] ssoName = ssoName[select] ssoKnown = ssoKnown[select] ssoCol = ssoCol[select,*] endif selectMOVIS=indgen(nbSSO) print, 'Number of asteroids in Bus-DeMeo:', nbSSO endelse ;--II.4--VISTA VHS Filters 7 Combinations ----------------------------------------------- filters=['Y', 'J', 'H','Ks'] nbFilt=n_elements(filters) nbComb = factorial(nbFilt)/ ( 2*factorial( nbFilt-2 )) label=['Y-J','Y-H','Y-Ks','J-H','J-Ks','H-Ks'] ; nbComb, 6 opciones <-traido de leeblos2 ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--- TAG --- III --- Compute Mean + StdDev for Each Class ------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--III.1--Output Structures ------------------------------------------------------------- covArr = replicate({correl:0., covar:0., angle:0., cos:0., sin:0., $ meanX:0., meanY:0., stdX:0., stdY:0., stdA:0., stdB:0.} ,nbClass,nbComb,nbComb) for kClass=0, nbClass-1 do begin ;--III.2--Select Only Asteroids from Current Class ----------------------------------- placeCl = where( strCmp(class, cLabel[kClass]), nbCl) print, cLabel[kClass]+'-types ('+strtrim(string(nbCl),2)+')' ;--III.3--Compute Correlations and Covariances --------------------------------------- locBDM = transpose(bdmCol[ placeCl,* ]) corr = correlate(locBDM) help, locBDM, corr cov = correlate(locBDM, /cov ) ;--III.3--Determine Could Orientation and Deviation ---------------------------------- for kComb1=0, nbComb-1 do begin ;--III.3.1--Sample Variance & Normalization of Covariance varComb1 = variance( bdmCol[placeCl,kComb1] ) cov[kComb1,*] /= varComb1 ;--III.3.2--Loop over Second Dimension for kComb2=kComb1+1, nbComb-1 do begin ;--III.3.3--Keep Correlation and Covariance covArr[kClass,kComb1,kComb2].correl = corr[kComb1,kComb2] covArr[kClass,kComb1,kComb2].covar = cov[kComb1,kComb2] ;--III.3.4--Compute Cloud Center covArr[kClass,kComb1,kComb2].meanX = mean( bdmCol[ placeCl,kComb1] ) covArr[kClass,kComb1,kComb2].meanY = mean( bdmCol[ placeCl,kComb2] ) covArr[kClass,kComb1,kComb2].stdX = stddev( bdmCol[ placeCl,kComb1] ) covArr[kClass,kComb1,kComb2].stdY = stddev( bdmCol[ placeCl,kComb2] ) ;--III.3.5--Compute Cloud Orientation covArr[kClass,kComb1,kComb2].angle = atan( covArr[kClass,kComb1,kComb2].covar ) covArr[kClass,kComb1,kComb2].cos = cos( covArr[kClass,kComb1,kComb2].angle ) covArr[kClass,kComb1,kComb2].sin = sin( covArr[kClass,kComb1,kComb2].angle ) ;--III.3.6--Compute Cloud Dispersion newX = bdmCol[placeCl,kComb1] - covArr[kClass,kComb1,kComb2].meanX newY = bdmCol[placeCl,kComb2] - covArr[kClass,kComb1,kComb2].meanY rotX = newX*covArr[kClass,kComb1,kComb2].cos + newY*covArr[kClass,kComb1,kComb2].sin rotY =-newX*covArr[kClass,kComb1,kComb2].sin + newY*covArr[kClass,kComb1,kComb2].cos covArr[kClass,kComb1,kComb2].stdA = stddev( rotX ) covArr[kClass,kComb1,kComb2].stdB = stddev( rotY ) endfor endfor endfor ;--III.3--Validation Figures ------------------------------------------------------------ if showGraphBDM ne 0 then begin ;--III.3.1--Loop over Filter Combinations for kComb1=0, nbComb-1 do begin for kComb2=kComb1+1, nbComb-1 do begin ;--III.3.2--Prepare Plot OutLay lonX = [ covArr[*,kComb1,kComb2].meanX+covArr[*,kComb1,kComb2].stdA, $ covArr[*,kComb1,kComb2].meanX-covArr[*,kComb1,kComb2].stdA, bdmCol[*,kComb1]] lonY = [ covArr[*,kComb1,kComb2].meanY+covArr[*,kComb1,kComb2].stdB, $ covArr[*,kComb1,kComb2].meanY-covArr[*,kComb1,kComb2].stdB, bdmCol[*,kComb2]] cgPlot, lonX, lonY, /NoData, xTitle=label[kComb1], yTitle=label[kComb2] ;--III.3.3--Plot Each Object within the Class for kClass=0, nbClass-1 do begin cur=where( strcmp(class,cLabel[kClass]) ) cgPlot, /OverPlot, bdmCol[cur,kComb1], bdmCol[cur,kComb2], $ psym=cSymb[kClass], color=cColor[kClass], symSize=sSymb[kClass] endfor ;--III.3.4--Plot Cross + Ellipse for each Class for kClass=0, nbClass-1 do begin cgPlot, covArr[kClass,kComb1,kComb2].meanX+covArr[kClass,kComb1,kComb2].stdA*cos(theta), $ covArr[kClass,kComb1,kComb2].meanY+covArr[kClass,kComb1,kComb2].stdB*sin(theta), color=cColor[kClass], /over cgErrPlot, covArr[kClass,kComb1,kComb2].meanX, color=cColor[kClass], $ covArr[kClass,kComb1,kComb2].meanY-covArr[kClass,kComb1,kComb2].stdB, $ covArr[kClass,kComb1,kComb2].meanY+covArr[kClass,kComb1,kComb2].stdB cgErrPlot, covArr[kClass,kComb1,kComb2].meanY, color=cColor[kClass], $ covArr[kClass,kComb1,kComb2].meanX-covArr[kClass,kComb1,kComb2].stdA, $ covArr[kClass,kComb1,kComb2].meanX+covArr[kClass,kComb1,kComb2].stdA, /horizontal cgText, covArr[kClass,kComb1,kComb2].meanX+covArr[kClass,kComb1,kComb2].stdA*0.3, $ covArr[kClass,kComb1,kComb2].meanY+covArr[kClass,kComb1,kComb2].stdB*0.3, cLabel[kClass], color=cColor[kClass] endfor stop endfor endfor endif ;--III.4--Validation Figures in a File -------------------------------------------------- if createGraphBDM eq 1 then begin ;--III.4--Validation Figures in a File -------------------------------------------------- ;--III.4.1-- Prepare plot cgPS_open, Filename=dirMOVIS+'DeMeo-Taxonomy.eps', /metric, /decomposed, /encapsulated, $ xSize=40, ySize=40, language_level=2, /quiet xLen = 0.176 yLen = xLen xSep = 0.01 ySep = 0.01 xBot = 0.085 yBot = 0.05 ;--III.4.2-- Loop over Combinations for kComb1=0, nbComb-1 do begin for kComb2=nbComb-1, kComb1+1, -1 do begin ;--III.4.3--Prepare Plot OutLay ;--III.4.3/A--XY Range lonX = [ covArr[*,kComb1,kComb2].meanX+covArr[*,kComb1,kComb2].stdA, $ covArr[*,kComb1,kComb2].meanX-covArr[*,kComb1,kComb2].stdA, bdmCol[*,kComb1]] lonY = [ covArr[*,kComb1,kComb2].meanY+covArr[*,kComb1,kComb2].stdB, $ covArr[*,kComb1,kComb2].meanY-covArr[*,kComb1,kComb2].stdB, bdmCol[*,kComb2]] xRang = ccRang[*,kComb1] yRang = ccRang[*,kComb2] ;--III.4.3/B--Graph position pos= [xBot + kComb1*xLen, $ yBot + (kComb2-1)*xLen, $ xBot + (kComb1+1)*xLen, $ yBot + (kComb2 )*xLen ] ;--III.4.3/C--XY Axes On the First column if kComb1 eq 0 then begin ;--Last line if kComb2 eq kComb1+1 then begin cgPlot, lonX, lonY, /NoData, /NoErase, /Normal, position=pos, $ xStyle=1, xRange=xRang, xTitle=label[kComb1], $ yStyle=1, yRange=yRang, yTitle=label[kComb2], $ xTickInt=ccTickInt[kComb1], xMinor=ccTickMin[kComb1], $ yTickInt=ccTickInt[kComb2], yMinor=ccTickMin[kComb2] ;--Other lines endif else begin cgPlot, lonX, lonY, /NoData, /NoErase, /Normal, position=pos, $ xStyle=1, xRange=xRang, xTickName = replicate(' ',20), $ yStyle=1, yRange=yRang, yTitle=label[kComb2], $ xTickInt=ccTickInt[kComb1], xMinor=ccTickMin[kComb1], $ yTickInt=ccTickInt[kComb2], yMinor=ccTickMin[kComb2] endelse ;--III.4.3/D--XY Axes On the other columns endif else begin ;--Last line if kComb2 eq kComb1+1 then begin cgPlot, lonX, lonY, /NoData, /NoErase, /Normal, position=pos, $ xStyle=1, xRange=xRang, xTitle=label[kComb1], $ yStyle=1, yRange=yRang, yTickName = replicate(' ',20), $ xTickInt=ccTickInt[kComb1], xMinor=ccTickMin[kComb1], $ yTickInt=ccTickInt[kComb2], yMinor=ccTickMin[kComb2] ;--Other lines endif else begin cgPlot, lonX, lonY, /NoData, /NoErase, /Normal, position=pos, $ xStyle=1, xRange=xRang, xTickName = replicate(' ',20), $ yStyle=1, yRange=yRang, yTickName = replicate(' ',20), $ xTickInt=ccTickInt[kComb1], xMinor=ccTickMin[kComb1], $ yTickInt=ccTickInt[kComb2], yMinor=ccTickMin[kComb2] endelse endelse ;--III.4.4--Plot Each Object within the Class for kClass=0, nbClass-1 do begin cur=where( strcmp(class,cLabel[kClass]) ) cgPlot, /OverPlot, bdmCol[cur,kComb1], bdmCol[cur,kComb2], $ psym=cSymb[kClass], color=cColor[kClass], symSize=sSymb[kClass] endfor ;--III.4.5--Plot Cross + Ellipse for each Class for kClass=0, nbClass-1 do begin ;--III.4.5/A--Long-Axis Line xA = [-1,1]*covArr[kClass,kComb1,kComb2].stdA yA = [0,0] rotXA = xA*covArr[kClass,kComb1,kComb2].cos - yA*covArr[kClass,kComb1,kComb2].sin + covArr[kClass,kComb1,kComb2].meanX rotYA = xA*covArr[kClass,kComb1,kComb2].sin + yA*covArr[kClass,kComb1,kComb2].cos + covArr[kClass,kComb1,kComb2].meanY cgPlot, rotXA, rotYA, color=cColor[kClass], /OverPlot ;--III.4.5/B--Short-Axis Line xB = [0,0] yB = [-1,1]*covArr[kClass,kComb1,kComb2].stdB rotXB = xB*covArr[kClass,kComb1,kComb2].cos - yB*covArr[kClass,kComb1,kComb2].sin + covArr[kClass,kComb1,kComb2].meanX rotYB = xB*covArr[kClass,kComb1,kComb2].sin + yB*covArr[kClass,kComb1,kComb2].cos + covArr[kClass,kComb1,kComb2].meanY cgPlot, rotXB, rotYB, color=cColor[kClass], /OverPlot ;--III.4.5/C--Ellipse x = covArr[kClass,kComb1,kComb2].stdA * cos(theta) y = covArr[kClass,kComb1,kComb2].stdB * sin(theta) rotX = x*covArr[kClass,kComb1,kComb2].cos - y*covArr[kClass,kComb1,kComb2].sin + covArr[kClass,kComb1,kComb2].meanX rotY = x*covArr[kClass,kComb1,kComb2].sin + y*covArr[kClass,kComb1,kComb2].cos + covArr[kClass,kComb1,kComb2].meanY cgPlot, rotX, rotY, color=cColor[kClass], /OverPlot endfor endfor endfor ;--III.4.5-- Display Key pos= [xBot + 3*xLen, $ yBot + (0)*yLen, $ xBot + (5)*xLen, $ yBot + (2 )*yLen ] cgPlot, [0,1], [0,1], /NoData, /NoErase, /Normal, position=pos, xStyle=5, yStyle=5 xKey = 0.2 xShi = 0.03 yKey = 0.90 yShi =-0.075 yTxt = -0.015 for kClass=0, nbClass-1 do begin cur=where( strcmp(class,cLabel[kClass]), nbCur ) cgPlot, xKey, yKey+kClass*yShi, psym=cSymb[kClass], color=cColor[kClass], symSize=sSymb[kClass], /OverPlot cgText, xKey+xShi, yKey+kClass*yShi+yTxt, cLabel[kClass]+'-types: ', color=cColor[kClass] cgText, xKey+xShi*13, yKey+kClass*yShi+yTxt, strtrim(string(nbCur),2), align=1, color=cColor[kClass] endfor cgPS_close, /png, Delete_PS=0 ;--III.5--Validation Figures in Many Files ---------------------------------------------- ;--III.5.1-- Loop over Combinations for kComb1=0, nbComb-1 do begin for kComb2=nbComb-1, kComb1+1, -1 do begin ;--III.5.2-- Prepare plot fileName = 'DeMeo-Taxonomy-'+label[kComb1]+'_vs_'+label[kComb2]+'.eps' cgPS_open, Filename=dirMOVIS+fileName, /metric, /decomposed, /encapsulated, $ xSize=15, ySize=15, language_level=2, /quiet ;--III.5.3--Prepare Plot OutLay xRang = ccRang[*,kComb1] yRang = ccRang[*,kComb2] cgPlot, 0,0, /NoData, charsize=1, $ xStyle=1, xRange=xRang, xTitle=label[kComb1], $ yStyle=1, yRange=yRang, yTitle=label[kComb2], $ xTickInt=ccTickInt[kComb1], xMinor=ccTickMin[kComb1], $ yTickInt=ccTickInt[kComb2], yMinor=ccTickMin[kComb2] ;--III.5.4--Plot Each Object within the Class for kClass=0, nbClass-1 do begin cur=where( strcmp(class,cLabel[kClass]) ) cgPlot, /OverPlot, bdmCol[cur,kComb1], bdmCol[cur,kComb2], $ psym=cSymb[kClass], color=cColor[kClass], symSize=sSymb[kClass] endfor ;--III.5.5--Plot Cross + Ellipse for each Class for kClass=0, nbClass-1 do begin ;--III.4.5/A--Long-Axis Line xA = [-1,1]*covArr[kClass,kComb1,kComb2].stdA yA = [0,0] rotXA = xA*covArr[kClass,kComb1,kComb2].cos - yA*covArr[kClass,kComb1,kComb2].sin + covArr[kClass,kComb1,kComb2].meanX rotYA = xA*covArr[kClass,kComb1,kComb2].sin + yA*covArr[kClass,kComb1,kComb2].cos + covArr[kClass,kComb1,kComb2].meanY cgPlot, rotXA, rotYA, color=cColor[kClass], /OverPlot ;--III.4.5/B--Short-Axis Line xB = [0,0] yB = [-1,1]*covArr[kClass,kComb1,kComb2].stdB rotXB = xB*covArr[kClass,kComb1,kComb2].cos - yB*covArr[kClass,kComb1,kComb2].sin + covArr[kClass,kComb1,kComb2].meanX rotYB = xB*covArr[kClass,kComb1,kComb2].sin + yB*covArr[kClass,kComb1,kComb2].cos + covArr[kClass,kComb1,kComb2].meanY cgPlot, rotXB, rotYB, color=cColor[kClass], /OverPlot ;--III.4.5/C--Ellipse x = covArr[kClass,kComb1,kComb2].stdA * cos(theta) y = covArr[kClass,kComb1,kComb2].stdB * sin(theta) rotX = x*covArr[kClass,kComb1,kComb2].cos - y*covArr[kClass,kComb1,kComb2].sin + covArr[kClass,kComb1,kComb2].meanX rotY = x*covArr[kClass,kComb1,kComb2].sin + y*covArr[kClass,kComb1,kComb2].cos + covArr[kClass,kComb1,kComb2].meanY cgPlot, rotX, rotY, color=cColor[kClass], /OverPlot cgText, covArr[kClass,kComb1,kComb2].meanX + 0.3*covArr[kClass,kComb1,kComb2].stdX, $ covArr[kClass,kComb1,kComb2].meanY + 0.3*covArr[kClass,kComb1,kComb2].stdY, cLabel[kClass], color=cColor[kClass] endfor cgPS_close, /png, Delete_PS=0 endfor endfor endif ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--- TAG --- IV --- Classify MOVIS Objects using BDM reference ------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; weightarr=fltarr(nbComb,nbComb)+1. ;--IV.1--Loop Over SSO in MOVIS --------------------------------------------------------- ssoClass = replicate('U',nbSSO) ssoComb = intarr(nbSSO) ;struc = fl ; struc = astelements(433) fname='~/W/Vista/datos/resumen.dat' ; to write down results from line 712 in file resumen.dat ;+++++++++++++++++++++ OPENW,1,fname ;+++++++++++++++++++++ printf, 1, 'Number', 'Name', 'Class', 'id','A','B','C','D','K','L','Q','S','T','X','V', $ 'Y-J','Y-H','Y-Ks','J-H ','J-Ks','H-Ks','a','e','i',$ format='(A-6,2x,A-20,6x,A-10,4x,11(A-2,6x),6(A-3,7x),3(A-1,10x))' ; struc = astelements(kSSO) for kSSO=0, nbSSO-1 do begin distArr=fltarr(nbClass) chi2Arr=fltarr(nbClass) testChiArr=fltarr(nbClass) ; Find DoF's +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ UsedFilter = intarr(4) ;***************************** >>>>>>>>>>>>>>>>>>>>>>> ; inside belows loop ;******************************************************* <<<<<<<<<<<<<<<<<<<<<<<< ;--IV.2--Loop Over Taxonomic Classes -------------------------------------------------- for kClass=0, nbClass-1 do begin nbUse = 0 chi2 = 0 ;--IV.3--Loop Over Combination of Filters ------------------------------------------- for kComb1=0, nbComb-1 do begin if ssoCol[kSSO,kComb1] ne -99.99 then $ for kComb2=round(nbComb)-1, kComb1+1, -1 do begin if ssoCol[kSSO,kComb2] ne -99.99 then begin ;;;;-------- aqui ; print , label[kcomb1], ' ', label[kcomb2] dX = ssoCol[kSSO,kComb1]-covArr[kClass,kComb1,kComb2].meanX dY = ssoCol[kSSO,kComb2]-covArr[kClass,kComb1,kComb2].meanY rotX = dX*covArr[kClass,kComb1,kComb2].cos + dY*covArr[kClass,kComb1,kComb2].sin rotY =-dX*covArr[kClass,kComb1,kComb2].sin + dY*covArr[kClass,kComb1,kComb2].cos chi2 += ( (rotX/covArr[kClass,kComb1,kComb2].stdA)^2. + $ (rotY/covArr[kClass,kComb1,kComb2].stdB)^2. ) * weightArr[kComb1,kComb2] nbUse++ ; print, ksso, kclass, label[kcomb1], label[kcomb2], dx, dy, $ ; format='(I5,2x,I2,2x,A-5,2x,A-5,3x,F8.3,2x,F8.3)' case label[kComb1] of 'Y-H' : fill = [0,2] ; 'H-Ks' : fill = [2,3] 'Y-Ks' : fill = [0,3] 'J-H' : fill = [1,2] 'J-Ks' : fill = [1,3] 'Y-J' : fill = [0,1] endcase UsedFilter[fill] = 1 ;print, UsedFilter case label[kComb2] of 'Y-H' : fill = [0,2] ; 'H-Ks' : fill = [2,3] 'Y-Ks' : fill = [0,3] 'J-H' : fill = [1,2] 'J-Ks' : fill = [1,3] 'Y-J' : fill = [0,1] endcase UsedFilter[fill] = 1 ;print, UsedFilter ;si tot filterused es 4 --> exit el loop if TOTAL(UsedFilter) eq 4 then begin kComb1 = nbComb kComb2 = nbComb endif endif endfor endfor DoF = TOTAL(UsedFilter)-1 ; if nBuse=0 means never used, look after ssocol[ksso,*] ne -99.99 ;this gives ne the correct filter pair ; idea: pair = where( ssocol[ksso,*] ne -99.99 ) ;nombre = label[pair] chi2Arr[kClass] = chi2 distArr[kClass] = sqrt(chi2) testChiArr[kClass] = 1 - chisqr_pdf( chi2, DoF ) ; dont forget to define it before endfor struc = astelements( ssoNum[kSSO] ) ; structure >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ;--IV.4--Find Most-Likely Class ------------------------------------------------------- ;yo ssoComb[kSSO] = nbUse if nbUse ge 1 then begin ;-TBD AQUI TIENES QUE TRABAJAR TU good = where( distArr le 5, nbClose ) ; antes era 5000, ahora 5 para ver los U mas cercanos, ; en promedio para todas las dimensiones de c/clase good = where( testChiArr ge 0.5, nbClose ) ; antes era 5000, ahora 5 para ver los U mas cercanos. Ultima version 50% if nbClose ge 1 then begin preOrder = sort( chi2Arr[good] ) order = good[preOrder] ; print, good ; print, preorder ; print, order ; stop ; used = where( ssoCol[kSSO,*] ne -99.99, nbUse ) ; ssoComb[kSSO] = nbUse ; ssoClass[kSSO] = cLabel[order[0]] ; dof = nbUse;-1 print, '' print, '['+strtrim(string(selectMOVIS[kSSO]+1),2)+']', ssoNum[kSSO], ssoName[kSSO], ssoKnown[kSSO], $ format='(A,4x,I6,2x,A-20,3x,"---[ ",A-1," ]---")' print, cLabel[order],format='(8x,11(A10))' print, 'dist :', distArr[order] , format='(A-6,2x,11(F10.3))' print, 'chi2 :', chi2Arr[order] , format='(A-6,2x,11(F10.3))' print, 'chi2R:', chi2Arr[order]/dof, format='(A-6,2x,11(F10.3))' ; chi2R : reduced chi2 print, 'DoF :', dof;, 2./dof print, 'testChi :', testChiArr[order] , format='(A-6,2x,11(F10.3))' printf, 1, ssoNum[kSSO], ssoName[kSSO], cLabel[order[0]], order[0], $ testChiArr[0], testChiArr[1], testChiArr[2], testChiArr[3], testChiArr[4],testChiArr[5], testChiArr[6],$ testChiArr[7], testChiArr[8], testChiArr[9], testChiArr[10], ssoCol[kSSO,0], ssoCol[kSSO,1], ssoCol[kSSO,2],$ ssoCol[kSSO,3], ssoCol[kSSO,4], ssoCol[kSSO,5], struc.a, struc.e,struc.i,$ ; ssoCol[kSSO,3], ssoCol[kSSO,4], ssoCol[kSSO,5], struc.a[kSSO], struc.e[kSSO],struc.i[kSSO],$ this only 1 ; format='(I6,2x,A-20,3x,"---[ ",A-1," ]---",4x,11(F6.4,2x),6(F8.4,2x),3(F8.4,2x))' format='(I6,", ",A-20,", ",A-1,", ",I2,", ",11(F6.4,", "),6(F8.4,", "),3(F8.4,", "))' ssoClass[kSSo] = cLabel[order[0]] ;; IDL> struc = astelements(433) ;; IDL> help, struc, /structure ;; ** Structure <1149d38>, 10 tags, length=72, data length=60, refs=1: ;; NUM LONG 433 ;; NAME STRING 'Eros' ;; H FLOAT 11.1600 ;; JD DOUBLE 2457600.5 ;; A DOUBLE 1.4579470 ;; E FLOAT 0.222606 ;; I FLOAT 10.8282 ;; O FLOAT 304.330 ;; W FLOAT 178.799 ;; N FLOAT 207.342 ;; ;; IDL> print, struc.num ;; 433 ;; IDL> print, struc.name ;; Eros ;; IDL> print, struc.a, struc.e, struc.i ;; 1.4579470 0.222606 10.8282 endif endif else begin ; endif else begin must be in the same line ;else begin para class 1 color *********** ;si if nBuse=0 means never used, look after ssocol[ksso,*] ne -99.99 ;this gives ne the correct filter pair pair = where( ssocol[ksso,*] ne -99.99 ) IDcol= label[pair] color = ssoCol[kSso,pair] case IDcol of 'Y-H': begin if color lt 0.45 then begin ssoClass[kSSO] = 'B' endif else begin if color gt 1.15 then ssoClass[kSSO] = 'A' ; if in one line doesn't req endif endelse end 'Y-J' : begin if color lt 0.25 then begin ssoClass[kSSO] = 'B' endif else begin if color gt 0.6 then ssoClass[kSSO] = 'V' endelse end 'Y-Ks' : begin if color lt 0.5 then begin ssoClass[kSSO] = 'B' endif else begin if color gt 1.3 then ssoClass[kSSO] = 'A' endelse end 'J-H' : begin if color lt 0.2 then begin ssoClass[kSSO] = 'V' endif else begin if color gt 0.65 then ssoClass[kSSO] = 'A' endelse end 'J-Ks' : begin if color lt 0.17 then begin ssoClass[kSSO] = 'V' endif else begin if color gt 0.8 then ssoClass[kSSO] = 'A' endelse end 'H-Ks' : begin if color gt -0.03 then begin ssoClass[kSSO] = 'V' endif else begin if color gt 0.2 then ssoClass[kSSO] = 'D' endelse end else: stop ;***************************;*************************** endcase endelse j2nextSSO: endfor CLOSE,1 ; ============================================================================== a ver ;stop cghistoplot, ssoComb, binsize=1, histdata=histComb, loca=xComb print, 'Number of asteroids in the catalog:', nbIn print, 'Number after uncertainty selection:', nbMOVIS forprint, xComb, histComb, textout=2 ; stop ;aqui where ssocomb igual 0 Y ssoclass no es 'U' sel0 = where( ssoComb eq 0 and ssoclass ne 'U', nbSel0 ) ;++++++++++++++++++++++++++++++ ;stop ; +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--- TAG --- V --- Selection of Valid Classification ------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--V.1--Only Select SSO with More than numSel colors ------------------------------------ sel = where( ssoComb ge numSel, nbSel ) print, 'Number with more than '+strtrim(string(numSel),2)+' colors:', nbSel if checkTaxo eq 0 then begin index = where( ssoComb ge numSel and $ ssoCol[*,0] ne -99.99 and $ ssoCol[*,1] ne -99.99 and $ ssoCol[*,2] ne -99.99 and $ ssoCol[*,3] ne -99.99 and $ ssoCol[*,4] ne -99.99 and $ ssoCol[*,5] ne -99.99 ) endif else begin index = indgen(nbSSO) endelse ;--V.2--Validation Figures -------------------------------------------------------------- if showGraphMOVIS eq 1 then begin ;--V.2.1--Loop over Filter Combinations for kComb1=0, nbComb-1 do begin for kComb2=kComb1+1, nbComb-1 do begin ;--V.2.2--Prepare Plot OutLay cgPlot, ssoCol[index,kComb1], ssoCol[index,kComb2], /NoData, xTitle=label[kComb1], yTitle=label[kComb2] ;--V.2.3--Plot Each Object within the Class for kClass=0, nbClass-1 do begin cur=where( strcmp(cLabel[kClass],ssoClass) ) cgPlot, /OverPlot, ssoCol[cur,kComb1], ssoCol[cur,kComb2], $ psym=cSymb[kClass], color=cColor[kClass], symSize=sSymb[kClass] endfor ;--V.2.4--Plot Cross + Ellipse for each Class for kClass=0, nbClass-1 do begin cgPlot, covArr[kClass,kComb1,kComb2].meanX+covArr[kClass,kComb1,kComb2].stdA*cos(theta), $ covArr[kClass,kComb1,kComb2].meanY+covArr[kClass,kComb1,kComb2].stdB*sin(theta), color=cColor[kClass], /over cgErrPlot, covArr[kClass,kComb1,kComb2].meanX, color=cColor[kClass], $ covArr[kClass,kComb1,kComb2].meanY-covArr[kClass,kComb1,kComb2].stdB, $ covArr[kClass,kComb1,kComb2].meanY+covArr[kClass,kComb1,kComb2].stdB cgErrPlot, covArr[kClass,kComb1,kComb2].meanY, color=cColor[kClass], $ covArr[kClass,kComb1,kComb2].meanX-covArr[kClass,kComb1,kComb2].stdA, $ covArr[kClass,kComb1,kComb2].meanX+covArr[kClass,kComb1,kComb2].stdA, /horizontal cgText, covArr[kClass,kComb1,kComb2].meanX+covArr[kClass,kComb1,kComb2].stdA*0.3, $ covArr[kClass,kComb1,kComb2].meanY+covArr[kClass,kComb1,kComb2].stdB*0.3, cLabel[kClass], color=cColor[kClass] endfor stop ; para probar los filtros que tengo endfor endfor endif ;--V.3--Validation Figures in a File -------------------------------------------------- if createGraphMOVIS eq 1 then begin ;--V.3.1-- Prepare plot if checkTaxo eq 0 then begin filename = 'MOVIS-Result-'+strtrim(string(numSel),2)+'-'+strtrim(string(medCut,format='(F4.2)'),2)+'.eps' endif else begin filename = 'DeMeo-Valid.eps' endelse ;medCut es el Median 0.12 en cette cas cgPS_open, Filename=dirMOVIS+filename, /metric, /decomposed, /encapsulated, $ xSize=40, ySize=40, language_level=2, /quiet xLen = 0.176 yLen = xLen xSep = 0.01 ySep = 0.01 xBot = 0.085 yBot = 0.05 ;--V.3.2-- Loop over Combinations for kComb1=0, nbComb-1 do begin for kComb2=nbComb-1, kComb1+1, -1 do begin ;--V.3.3--Prepare Plot OutLay ;--V.3.3/A--XY Range lonX = [ covArr[*,kComb1,kComb2].meanX+covArr[*,kComb1,kComb2].stdA, $ covArr[*,kComb1,kComb2].meanX-covArr[*,kComb1,kComb2].stdA, ssoCol[index,kComb1]] lonY = [ covArr[*,kComb1,kComb2].meanY+covArr[*,kComb1,kComb2].stdB, $ covArr[*,kComb1,kComb2].meanY-covArr[*,kComb1,kComb2].stdB, ssoCol[index,kComb2]] xRang = ccRang[*,kComb1] yRang = ccRang[*,kComb2] ;--V.3.3/B--Graph position pos= [xBot + kComb1*xLen, $ yBot + (kComb2-1)*xLen, $ xBot + (kComb1+1)*xLen, $ yBot + (kComb2 )*xLen ] ;--V.3.3/C--XY Axes On the First column if kComb1 eq 0 then begin ;--Last line if kComb2 eq kComb1+1 then begin cgPlot, lonX, lonY, /NoData, /NoErase, /Normal, position=pos, $ xStyle=1, xRange=xRang, xTitle=label[kComb1], $ yStyle=1, yRange=yRang, yTitle=label[kComb2] ;--Other lines endif else begin cgPlot, lonX, lonY, /NoData, /NoErase, /Normal, position=pos, $ xStyle=1, xRange=xRang, xTickName = replicate(' ',20), $ yStyle=1, yRange=yRang, yTitle=label[kComb2] endelse ;--V.3.3/D--XY Axes On the other columns endif else begin ;--Last line if kComb2 eq kComb1+1 then begin cgPlot, lonX, lonY, /NoData, /NoErase, /Normal, position=pos, $ xStyle=1, xRange=xRang, xTitle=label[kComb1], $ yStyle=1, yRange=yRang, yTickName = replicate(' ',20) ;--Other lines endif else begin cgPlot, lonX, lonY, /NoData, /NoErase, /Normal, position=pos, $ xStyle=1, xRange=xRang, xTickName = replicate(' ',20), $ yStyle=1, yRange=yRang, yTickName = replicate(' ',20) endelse endelse ;-u cur=where( strcmp('U',ssoClass), nbCur ) ; ++++++++++++++++++++++++= cgPlot, /OverPlot, ssoCol[cur,kComb1], ssoCol[cur,kComb2], $ psym=cSymb[nbClass], color=cColor[nbClass], symSize=sSymb[nbClass] ;++++++++++++++++++ def symb col size print, 'HACIENDO LOS U!!!!! ' + filename print, nbcur ;--V.3.4--Plot Each Object within the Class for kClass=0, nbClass-1 do begin cur=where( strcmp(cLabel[kClass],ssoClass) and ssoComb ge numSel, nbCur ) ; print, label[kComb1], label[kComb2], cLabel[kClass], nbCur, format='(A-5,A-5,A-5,I5)' cgPlot, /OverPlot, ssoCol[cur,kComb1], ssoCol[cur,kComb2], $ psym=cSymb[kClass], color=cColor[kClass], symSize=sSymb[kClass] endfor ;--V.3.5--Plot Cross + Ellipse for each Class for kClass=0, nbClass-1 do begin ;--V.3.5/A--Long-Axis Line xA = [-1,1]*covArr[kClass,kComb1,kComb2].stdA yA = [0,0] rotXA = xA*covArr[kClass,kComb1,kComb2].cos - yA*covArr[kClass,kComb1,kComb2].sin + covArr[kClass,kComb1,kComb2].meanX rotYA = xA*covArr[kClass,kComb1,kComb2].sin + yA*covArr[kClass,kComb1,kComb2].cos + covArr[kClass,kComb1,kComb2].meanY cgPlot, rotXA, rotYA, color=cColor[kClass], /OverPlot ;--V.3.5/B--Short-Axis Line xB = [0,0] yB = [-1,1]*covArr[kClass,kComb1,kComb2].stdB rotXB = xB*covArr[kClass,kComb1,kComb2].cos - yB*covArr[kClass,kComb1,kComb2].sin + covArr[kClass,kComb1,kComb2].meanX rotYB = xB*covArr[kClass,kComb1,kComb2].sin + yB*covArr[kClass,kComb1,kComb2].cos + covArr[kClass,kComb1,kComb2].meanY cgPlot, rotXB, rotYB, color=cColor[kClass], /OverPlot ;--V.3.5/C--Ellipse x = covArr[kClass,kComb1,kComb2].stdA * cos(theta) y = covArr[kClass,kComb1,kComb2].stdB * sin(theta) rotX = x*covArr[kClass,kComb1,kComb2].cos - y*covArr[kClass,kComb1,kComb2].sin + covArr[kClass,kComb1,kComb2].meanX rotY = x*covArr[kClass,kComb1,kComb2].sin + y*covArr[kClass,kComb1,kComb2].cos + covArr[kClass,kComb1,kComb2].meanY cgPlot, rotX, rotY, color=cColor[kClass], /OverPlot endfor endfor endfor ;--V.3.6-- Display Key pos= [xBot + 3*xLen, $ yBot + (0)*yLen, $ xBot + (5)*xLen, $ yBot + (2 )*yLen ] cgPlot, [0,1], [0,1], /NoData, /NoErase, /Normal, position=pos, xStyle=5, yStyle=5 xKey = 0.2 xShi = 0.03 yKey = 0.90 yShi =-0.075 yTxt = -0.015 for kClass=0, nbClass-1 do begin cur=where( strcmp(cLabel[kClass],ssoClass) and ssoComb ge numSel, nbCur ) ; cur=where( strcmp(ssoClass,cLabel[kClass]), nbCur ) cgPlot, xKey, yKey+kClass*yShi, psym=cSymb[kClass], color=cColor[kClass], symSize=sSymb[kClass], /OverPlot cgText, xKey+xShi, yKey+kClass*yShi+yTxt, cLabel[kClass]+'-types: ', color=cColor[kClass] cgText, xKey+xShi*13, yKey+kClass*yShi+yTxt, strtrim(string(nbCur),2), align=1, color=cColor[kClass] endfor cgPS_close, /png, Delete_PS=0 ;--V.4--Validation Figures in Many Files ---------------------------------------------- ;--V.4.1-- Loop over Combinations for kComb1=0, nbComb-1 do begin for kComb2=nbComb-1, kComb1+1, -1 do begin ;--V.4.2-- Prepare plot if checkTaxo eq 0 then begin filename = 'MOVIS-Result-'+strtrim(string(numSel),2)+'-'+strtrim(string(medCut,format='(F4.2)'),2)+$ '-'+label[kComb1]+'_vs_'+label[kComb2]+'.eps' endif else begin fileName = 'DeMeo-Valid-'+label[kComb1]+'_vs_'+label[kComb2]+'.eps' endelse cgPS_open, Filename=dirMOVIS+fileName, /metric, /decomposed, /encapsulated, $ xSize=15, ySize=15, language_level=2, /quiet ;--V.4.3--Prepare Plot OutLay xRang = ccRang[*,kComb1] yRang = ccRang[*,kComb2] cgPlot, 0,0, /NoData, charsize=1, $ xStyle=1, xRange=xRang, xTitle=label[kComb1], $ yStyle=1, yRange=yRang, yTitle=label[kComb2], $ xTickInt=ccTickInt[kComb1], xMinor=ccTickMin[kComb1], $ yTickInt=ccTickInt[kComb2], yMinor=ccTickMin[kComb2] ; cur=where( strcmp('U',ssoClass), nbCur ) ; ++++++++++++++++++++++++= cgPlot, /OverPlot, ssoCol[cur,kComb1], ssoCol[cur,kComb2], $ psym=cSymb[nbClass], color=cColor[nbClass], symSize=sSymb[nbClass] ;++++++++++++++++++ def symb col size ;--V.4.4--Plot Each Object within the Class for kClass=0, nbClass-1 do begin cur=where( strcmp(cLabel[kClass],ssoClass) and ssoComb ge numSel, nbCur ) ; cur=where( strcmp(cLabel[kClass],ssoClass), nbCur ) ; print, label[kComb1], label[kComb2], cLabel[kClass], nbCur, format='(A-5,A-5,A-5,I5)' cgPlot, /OverPlot, ssoCol[cur,kComb1], ssoCol[cur,kComb2], $ psym=cSymb[kClass], color=cColor[kClass], symSize=sSymb[kClass] endfor ;--V.4.5--Plot Cross + Ellipse for each Class for kClass=0, nbClass-1 do begin ;--V.4.5/A--Long-Axis Line xA = [-1,1]*covArr[kClass,kComb1,kComb2].stdA yA = [0,0] rotXA = xA*covArr[kClass,kComb1,kComb2].cos - yA*covArr[kClass,kComb1,kComb2].sin + covArr[kClass,kComb1,kComb2].meanX rotYA = xA*covArr[kClass,kComb1,kComb2].sin + yA*covArr[kClass,kComb1,kComb2].cos + covArr[kClass,kComb1,kComb2].meanY cgPlot, rotXA, rotYA, color=cColor[kClass], /OverPlot ;--V.4.5/B--Short-Axis Line xB = [0,0] yB = [-1,1]*covArr[kClass,kComb1,kComb2].stdB rotXB = xB*covArr[kClass,kComb1,kComb2].cos - yB*covArr[kClass,kComb1,kComb2].sin + covArr[kClass,kComb1,kComb2].meanX rotYB = xB*covArr[kClass,kComb1,kComb2].sin + yB*covArr[kClass,kComb1,kComb2].cos + covArr[kClass,kComb1,kComb2].meanY cgPlot, rotXB, rotYB, color=cColor[kClass], /OverPlot ;--V.4.5/C--Ellipse x = covArr[kClass,kComb1,kComb2].stdA * cos(theta) y = covArr[kClass,kComb1,kComb2].stdB * sin(theta) rotX = x*covArr[kClass,kComb1,kComb2].cos - y*covArr[kClass,kComb1,kComb2].sin + covArr[kClass,kComb1,kComb2].meanX rotY = x*covArr[kClass,kComb1,kComb2].sin + y*covArr[kClass,kComb1,kComb2].cos + covArr[kClass,kComb1,kComb2].meanY cgPlot, rotX, rotY, color=cColor[kClass], /OverPlot cgText, covArr[kClass,kComb1,kComb2].meanX + 0.3*covArr[kClass,kComb1,kComb2].stdX, $ covArr[kClass,kComb1,kComb2].meanY + 0.3*covArr[kClass,kComb1,kComb2].stdY, cLabel[kClass], color=cColor[kClass] endfor cgPS_close, /png, Delete_PS=0 endfor endfor endif ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--- TAG --- VI --- Classification Validation with DeMeo ------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; j2xm: if checkTaxo eq 1 then begin ;ben do comments & all comArr = intarr(nbClass,nbClass)+2 comArr[indgen(nbClass),indgen(nbClass)] = 0 ok_A=[6,7,10] & A=0 & comArr[A,ok_A]=1 & comArr[ok_A,A]=1 ok_B=[2,9] & B=1 & comArr[B,ok_B]=1 & comArr[ok_B,B]=1 ok_C=[9] & C=2 & comArr[C,ok_C]=1 & comArr[ok_C,C]=1 ok_D=[8] & D=3 & comArr[D,ok_D]=1 & comArr[ok_D,D]=1 ok_K=[5,7] & K=4 & comArr[K,ok_K]=1 & comArr[ok_K,K]=1 ok_L=[7] & L=5 & comArr[L,ok_L]=1 & comArr[ok_L,L]=1 ok_Q=[7,10] & Q=6 & comArr[Q,ok_Q]=1 & comArr[ok_Q,Q]=1 ok_S=[10] & S=7 & comArr[S,ok_S]=1 & comArr[ok_S,S]=1 ; ok_T=[8] & T=8 & comArr[T,ok_T]=1 & comArr[ok_T,T]=1 ok_X=[2] & X=9 & comArr[X,ok_X]=1 & comArr[ok_X,X]=1 ok_V=[7] & V=10 & comArr[V,ok_V]=1 & comArr[ok_V,V]=1 xmN = intarr(nbClass,nbClass) xmF = fltarr(nbClass,nbClass) xmL = strarr(nbClass,nbClass) xmV = fltarr(nbClass,nbClass) for kClass=0, nbClass-1 do begin cur = where( strcmp( Class, cLabel[kClass]), nbValid ) local = ssoClass[cur] uniqClass = local[ uniq( local, sort(local) ) ] nbUniq=n_elementS(uniqClass) numPerClass=intarr(nbUniq) frqPerClass=fltarr(nbUniq) for kU=0, nbUniq-1 do begin pp = where( strcmp(local,uniqClass[kU]), nbPP ) numPerClass[kU]=nbPP endfor frqPerClass=100.*numPerClass/nbValid ord=reverse(sort(numPerClass)) if strcmp( cLabel[kClass], 'T' ) then begin ord=sort(numPerClass) endif xmN[kClass,0:nbUniq-1] = numPerClass[ord] xmF[kClass,0:nbUniq-1] = frqPerClass[ord] xmL[kClass,0:nbUniq-1] = uniqClass[ord] print, '' print, cLabel[kClass], ' -- ', uniqClass[ord[0]], numPerClass[ord[0]], frqPerClass[ord[0]], $ format='(A1,A4,A1,3x,I3,3x,F6.2)' if nbUniq gt 1 then begin for kU=1, nbUniq-1 do begin print, ' ', ' ', uniqClass[ord[kU]], numPerClass[ord[kU]], frqPerClass[ord[kU]], $ format='(A1,A4,A1,3x,I3,3x,F6.2)' endfor print, ' ', ' ', 'Tot', nbValid, total(frqPerClass), $ format='(A1,A3,A3,2x,I3,3x,F6.2)' endif endfor statValid=fltarr(nbClass,3) statValid[*,0] = xmF[*,0] cgPS_open, Filename=dirMOVIS+'DeMeo-Valid-Fraction.eps', /metric, /decomposed, /encapsulated, $ xSize=25, ySize=18, language_level=2, /quiet cgBarPlot, reverse(xmF[*,0]), /Rotate, /normal, $ xStyle=1, xRange=[0,100], $ yStyle=1, yRange=[0,nbClass-1], $ barNames=reverse(cLabel), colors=reverse(cColor), $ xTitle='Fraction of classification', $ yTitle='VISTA classification', $ barCoords=yArr, $ position=[0.08,0.12,0.97,0.97] barSize=0.6675 barOffset=-0.33 last = xmF[*,0] for kE=1, nbClass-1 do begin colorVec=replicate('Black',nbClass) ; for kClass=nbClass-1,0,-1 do begin for kClass=0, nbClass-1 do begin p=where( strcmp(cLabel, xmL[kClass,kE]) ) if p[0] ne -1 then begin colorVec[kClass] = cColor[p[0]] if comArr[kClass,p[0]] lt 2 then xmV[kClass,kE]=1 statValid[*, comArr[kClass,p[0]]] += xmF[*,kClass] if xmV[kClass,kE] eq 1 and xmF[kClass,kE] gt 0.5 then begin base = last last[kClass] += xmF[kClass,kE] xx = base[kClass] + [0,xmF[kClass,kE],xmF[kClass,kE],0,0] yy = [0,0,1,1,0]*barSize + barOffSet + yArr[nbClass-1-kClass] cgPlot, /over, xx, yy, color=colorVec[kClass] cgText, mean(xx[0:1]), yy[1] + (yy[2]-yy[1])/4., align=0.5, xmL[kClass,kE], color=colorVec[kClass] endif endif endfor endfor cgPS_close, /png, Delete_PS=0 endif ; bad = where( strcmp(ssoClass ,'U'),nbu ) ; print, nbu, nbsso, nbSSO-nbU ; forprint, movis.num, movis.name, ssoClass, ssoComb, textout=2, $ ; format='(I6,1x,A-20,1x,A1,1x,I1)', subset=sel ; bad = where( strcmp(ssoClass ,'U'),nbu ) ; printf, nbu, nbsso, nbSSO-nbU ; # , name, class, [ 11 x probCi ] , Y-J , Y-H ; fname='~/W/Vista/datos/resumen.dat' ; ; OPENW,1,fname ; ; ; printf, '' ; printf,1, '['+strtrim(string(selectMOVIS[kSSO]+1),2)+']', ssoNum[kSSO], ssoName[kSSO], ssoKnown[kSSO] ;, testChiArr[kSSO] $ ; format='(A,4x,I6,2x,A-20,3x,"---[ ",A-3," ]---")' ; ; CLOSE,1 ; printf, 'testChi :', testChiArr[order] ; A B C D K L Q S ; T X V ; apres ajouter a, e , i end