;-como fnciona where a=indgen(5) ;-como fnciona where select = where( (a mod 2) eq 0, nbSel ) ;-como fnciona where ;-como fnciona where print, a ;-como fnciona where ;-como fnciona where for kSel=0, nbSel-1 do begin ;-como fnciona where ;-como fnciona where print, kSel, select[kSel], a[select[kSel]] ;-como fnciona where ;-como fnciona where ;-como fnciona where endfor ;-como fnciona where ;-como fnciona where ;-como fnciona where; print, select ;-como fnciona where; print, '' ;-como fnciona where print, a[select] ;-como fnciona where ;-como fnciona where ;-como fnciona wherestop ;-First Guess errCutYJ = 0.25 errCutYH = 0.35 errCutYK = 0.35 errCutJH = 0.25 errCutJK = 0.25 errCutHK = 0.15 ;-More stringent errCutYJ = 0.1 errCutYH = 0.1 errCutYK = 0.1 errCutJH = 0.1 errCutJK = 0.1 errCutHK = 0.1 medCut = median( [errCutYJ,errCutYH,errCutYK,errCutJH,errCutJK,errCutHK] ) ;-What to do showGraphBDM = 0 showGraphMOVIS = 0 checkTaxo = 1 ;-User choice numSel = 6 dirMOVIS = '/home/bcarry/data/mining/movis/' ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--- TAG --- I --- Define classes ------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--I.1--Class Labels -------------------------------------------------------------------- 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 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 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 ;--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' ) ;--II.2.2--Selection based on threshold select = 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, nbSel) ;--II.2.3--Selection based on threshold movis = movis[select] 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 endif else begin ssoCol = bdmCol nbSSO = n_elements( class ) 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 ------------------------------------------------------------- meanArr = fltarr(nbClass,nbComb) stdArr = fltarr(nbClass,nbComb) ;--III.2--Loop over Classes ------------------------------------------------------------- for kClass=0, nbClass-1 do begin ;--III.2.1--Search Asteroids of current Class placeCl = where( strCmp(class, cLabel[kClass]), nbCl) print, cLabel[kClass]+'-types ('+strtrim(string(nbCl),2)+')' ;--III.2.2--Loop over Filter Combination ;-TBD: compute correlation & orientation for kComb=0, nbComb-1 do begin meanArr[kClass,kComb] = mean( bdmCol[ placeCl,kComb] ) stdArr[ kClass,kComb] = stddev( bdmCol[ placeCl,kComb] ) endfor endfor ;--III.3--Validation Figures on Screen --------------------------------------------------- if showGraphBDM eq 1 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 = [ meanArr[*,kComb1]+stdArr[*,kComb1], meanArr[*,kComb1]-stdArr[*,kComb1], bdmCol[*,kComb1]] lonY = [ meanArr[*,kComb2]+stdArr[*,kComb2], meanArr[*,kComb2]-stdArr[*,kComb2], 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, meanArr[kClass,kComb1]+stdArr[kClass,kComb1]*cos(theta), $ meanArr[kClass,kComb2]+stdArr[kClass,kComb2]*sin(theta), color=cColor[kClass], /over cgErrPlot, meanArr[kClass,kComb1], color=cColor[kClass], $ meanArr[kClass,kComb2]-stdArr[kClass,kComb2], meanArr[kClass,kComb2]+stdArr[kClass,kComb2] cgErrPlot, meanArr[kClass,kComb2], color=cColor[kClass], meanArr[kClass,kComb1]-stdArr[kClass,kComb1], $ meanArr[kClass,kComb1]+stdArr[kClass,kComb1], /horizontal cgText, meanArr[kClass,kComb1]+stdArr[kClass,kComb1]*0.3, $ meanArr[kClass,kComb2]+stdArr[kClass,kComb2]*0.3, cLabel[kClass], color=cColor[kClass] endfor stop endfor endfor endif ;--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 = [ meanArr[*,kComb1]+stdArr[*,kComb1], meanArr[*,kComb1]-stdArr[*,kComb1], bdmCol[*,kComb1] ] lonY = [ meanArr[*,kComb2]+stdArr[*,kComb2], meanArr[*,kComb2]-stdArr[*,kComb2], bdmCol[*,kComb2] ] xRang = [min(lonX), max(lonX)] yRang = [min(lonY), max(lonY)] ;--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] ;--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 ;--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) ;--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 ;--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 cgPlot, meanArr[kClass,kComb1]+stdArr[kClass,kComb1]*cos(theta), $ meanArr[kClass,kComb2]+stdArr[kClass,kComb2]*sin(theta), color=cColor[kClass], /over cgErrPlot, meanArr[kClass,kComb1], color=cColor[kClass], $ meanArr[kClass,kComb2]-stdArr[kClass,kComb2], meanArr[kClass,kComb2]+stdArr[kClass,kComb2] cgErrPlot, meanArr[kClass,kComb2], color=cColor[kClass], meanArr[kClass,kComb1]-stdArr[kClass,kComb1], $ meanArr[kClass,kComb1]+stdArr[kClass,kComb1], /horizontal cgText, meanArr[kClass,kComb1]+stdArr[kClass,kComb1]*0.3, $ meanArr[kClass,kComb2]+stdArr[kClass,kComb2]*0.3, cLabel[kClass], color=cColor[kClass] endfor endfor endfor cgPS_close, /png, Delete_PS=0 ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--- TAG --- IV --- Classify MOVIS Objects using BDM reference ------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--IV.1--Loop Over SSO in MOVIS --------------------------------------------------------- ssoClass = strarr(nbSSO) ssoComb = intarr(nbSSO) for kSSO=0, nbSSO-1 do begin distArr=fltarr(nbClass) chi2Arr=fltarr(nbClass) ;--IV.2--Loop Over Taxonomic Classes -------------------------------------------------- for kClass=0, nbClass-1 do begin distance = 0. chi2=0. nbUse = 0 ;--IV.3--Loop Over Combination of Filters ------------------------------------------- for kComb=0, nbComb-1 do begin ; sobre todas las combinaciones... nComb=6 ;--IV.3.1-- Average and StdDev of current Class c = meanArr[kClass,kComb] sc = stdArr[kClass,kComb] ;--IV.3.2-- IF SSO has color in current Combination -> add distance if ssoCol[kSSO,kComb] ne -99.99 then begin chi2 += ( (ssoCol[kSSO,kComb]-c)/sc )^2 distance += ( (ssoCol[kSSO,kComb]-c)/sc )^2 nbUse++ endif endfor ;--IV.3--Loop Over Combination of Filters ------------------------------------------- rms = sqrt( distance/nbUse ) distArr[kClass] = rms chi2Arr[kClass] = chi2 endfor ;--IV.4--Find Most-Likely Class ------------------------------------------------------- order = sort( distArr ) ssoComb[kSSO] = nbUse ssoClass[kSSO] = cLabel[order[0]] ;TBD: Check different of distARR --> select ALL classes with similar distance ;- give preference to main classes over minor classes order = sort( chi2Arr ) ssoComb[kSSO] = nbUse ssoClass[kSSO] = cLabel[order[0]] print, '' print, kSSO, '['+class[kSSO]+']', cLabel[order] print, chi2Arr[order] dof = nbUse-1 print, dof cut90 = [0 , 2.71, 4.60, 6.25, 7.78, 9.24] cut95 = [0 , 3.84, 5.88, 7.82, 9.49, 11.07] cut99 = [0 , 6.64, 9.21, 11.34, 13.28, 15.09] sel90 = where( chi2Arr[order] le cut90[dof], nb90 ) sel95 = where( chi2Arr[order] le cut95[dof], nb95 ) sel99 = where( chi2Arr[order] le cut99[dof], nb99 ) if nb90 gt 0 then begin print, 'Sel 90: ', cLabel[order[sel90]] print, ' ', chi2Arr[order[sel90]] endif ; if nb95 gt 0 then begin ; print, 'Sel 95: ', cLabel[order[sel95]] ; print, ' ', chi2Arr[order[sel95]] ; endif ; ; if nb99 gt 0 then begin ; print, 'Sel 99: ', cLabel[order[sel99]] ; print, ' ', chi2Arr[order[sel99]] ; endif ;stop endfor ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--- TAG --- V --- Selection of Valid Classification ------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--V.1--Only Select SSO with More than numSel colors ------------------------------------ sel = where( ssoComb ge numSel, nbSel ) print, 'Number of asteroids 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, meanArr[kClass,kComb1]+stdArr[kClass,kComb1]*cos(theta), $ meanArr[kClass,kComb2]+stdArr[kClass,kComb2]*sin(theta), color=cColor[kClass], /over cgErrPlot, meanArr[kClass,kComb1], color=cColor[kClass], $ meanArr[kClass,kComb2]-stdArr[kClass,kComb2], meanArr[kClass,kComb2]+stdArr[kClass,kComb2] cgErrPlot, meanArr[kClass,kComb2], color=cColor[kClass], meanArr[kClass,kComb1]-stdArr[kClass,kComb1], $ meanArr[kClass,kComb1]+stdArr[kClass,kComb1], /horizontal cgText, meanArr[kClass,kComb1]+stdArr[kClass,kComb1]*0.3, $ meanArr[kClass,kComb2]+stdArr[kClass,kComb2]*0.3, cLabel[kClass], color=cColor[kClass] endfor stop endfor endfor endif ;--V.3--Validation Figures in a File -------------------------------------------------- ;--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 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 = [ meanArr[*,kComb1]+stdArr[*,kComb1], meanArr[*,kComb1]-stdArr[*,kComb1], ssoCol[index,kComb1]] lonY = [ meanArr[*,kComb2]+stdArr[*,kComb2], meanArr[*,kComb2]-stdArr[*,kComb2], ssoCol[index,kComb2]] xRang = [min(lonX), max(lonX)] yRang = [min(lonY), max(lonY)] ;--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 ;--V.3.4--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.3.5--Plot Cross + Ellipse for each Class for kClass=0, nbClass-1 do begin cgPlot, meanArr[kClass,kComb1]+stdArr[kClass,kComb1]*cos(theta), $ meanArr[kClass,kComb2]+stdArr[kClass,kComb2]*sin(theta), color=cColor[kClass], /over cgErrPlot, meanArr[kClass,kComb1], color=cColor[kClass], $ meanArr[kClass,kComb2]-stdArr[kClass,kComb2], meanArr[kClass,kComb2]+stdArr[kClass,kComb2] cgErrPlot, meanArr[kClass,kComb2], color=cColor[kClass], meanArr[kClass,kComb1]-stdArr[kClass,kComb1], $ meanArr[kClass,kComb1]+stdArr[kClass,kComb1], /horizontal cgText, meanArr[kClass,kComb1]+stdArr[kClass,kComb1]*0.3, $ meanArr[kClass,kComb2]+stdArr[kClass,kComb2]*0.3, cLabel[kClass], color=cColor[kClass] endfor endfor endfor cgPS_close, /png, Delete_PS=0 ; forprint, movis.num, movis.name, ssoClass, ssoComb, textout=2, $ ; format='(I6,1x,A-20,1x,A1,1x,I1)', subset=sel end