onlyValidFig=0 verbose = 1 ;-Processing: classTaxo = 0 ;-To class Bus-DeMeo sample findAEI = 01 ;-Make it MUCH slower ;-figures: plotReferenceMulti = 0 ;-Make a plot of Reference Sample - all Filters in a single File plotReferenceSingle = 0 ;-Make a plot of Reference Sample - Each Filter in a single File ;--- probability selection minProba = 0.5 ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--- TAG --- I --- General User Choices ------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; if verbose eq 1 then begin print, '' print, '--------------------------------------------------------------------------------' print, ' I. Load Configuration' endif ;--I.1-- Bunch of Filenames -------------------------------------------------------------; fileColorCut = 'colors.cut' fileClassDef = 'class.def' ;--I.2-- Directory ----------------------------------------------------------------------; ;--I.2.1-- Root Directory spawn, 'hostname', host case strTrim(host) of 'hyperion': root = '/data/colors/movis/' 'endymion': root = '/home/bcarry/data/mining/movis/' 'gnanzi': root = '/home/mpajuelo/W/Vista/' 'IMCCE5': root = 'PON_TU_PATH' else: stop endcase ;--I.2.2-- SubTree dirDATA = root + 'data/' dirFIG = root + 'figures/' dirTAXO = root + 'taxonomy/' dirPAR = root + 'param/' ;--I.2.3-- Create Directories if Needed if not file_test(dirDATA,/Directory) then file_mkdir, dirDATA if not file_test(dirFIG ,/Directory) then file_mkdir, dirFIG if not file_test(dirTAXO,/Directory) then file_mkdir, dirTAXO if not file_test(dirPAR ,/Directory) then file_mkdir, dirPAR ;--I.3-- VISTA-Related Parameters -------------------------------------------------------; colorCut = movisReadParam( dirPAR+fileColorCut, /ColorCut ) medCut = median( [colorCut.YJ,colorCut.YH,colorCut.YK,colorCut.JH,colorCut.JK,colorCut.HK] ) ;--I.4-- Asteroid Taxonomy --------------------------------------------------------------; ;--I.4.1-- Read Classes Colors, Symbols, ... classDef = movisReadParam( dirPAR+fileClassDef, /Class ) validClass = where( classDef.use eq 1, nbClass ) class = classDef[validClass] ;--I.4.2-- GUI: Classes print, 'List of classes considered here ('+strtrim(string(nbClass),2)+') ' print, class.label ;--I.4.3-- Minimum number of colors to include minNumOfColor = 1 ;--I.5-- VISTA VHS Combination of Filters -----------------------------------------------; filters=['Y', 'J', 'H','Ks'] nbFilt=n_elements(filters) nbComb = factorial(nbFilt)/ ( 2*factorial( nbFilt-2 )) combLabel=['Y-J','Y-H','Y-Ks','J-H','J-Ks','H-Ks'] ;--I.x-- to be changed to config files --------------------------------------------------; ;-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 theta=2*!PI*findgen(361)/360 ;--- class colors 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 ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--- TAG --- II --- Read Learning Data Sets ------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; if verbose eq 1 then begin print, '' print, '--------------------------------------------------------------------------------' print, ' II. Read Learning Data Sets' endif ;--II.1-- Bus-DeMeo Spectra Converted to 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, dirDATA+'bdm-col.csv',f=data,$ bdmClass,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]] ;-TBD - Redo Bus-DeMeo to VHS filters, with MOVIS-C format!!!!!! ; 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-- Extract MOVIS Objects with Known Taxonomy -------------------------------------; ;-TBD - 1 read movis ;-TBD - 2 read known classes ;-TBD - 3 crossmatch ;-TBD - 4 fill simple array with colors ; movis = vistaReadMOVIS( dirDATA+'MOVIS-C.csv' ) ;--II.3-- SMASS Asteroids from Binzel ---------------------------------------------------; ;-TBD - 1 read list of known taxo class ;-TBD - 2 read corresponding asteroids ;-TBD - 3 convert spectra to colors ;-TBD - 4 write refernece file ;-TBD ---------------- Use convertTaxoToColor.pro !!!!!!!! ;--II.x-- Concatenate Learning Samples --------------------------------------------------; ;-TBD - add keyword for using BDM, SMASS, MOVIS -> concatenate ;-TBD - Use num:name: tags in ref structure, see sso in IV below nbRef = n_elements(bdmClass) emptySSO={num:0L, name:'', class:'', $ ;-Basic Id: Number, Name, Taxonomic Class N:0, col:fltarr(6), id:-1, $ ;-Obs.: Number of filter, colors, -CLASS ID- a:0.d, e:0., i:0., H:0., $ ;-Dyn.: Orbital elements proba:fltarr(nbClass) } ;-Taxo: Probability for each Class ref = replicate(emptySSO,nbRef) ref.class= bdmClass for k=0,nbRef-1 do ref[k].col = bdmCol[k,*] ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--- TAG --- III --- Define Class Ellipses (mean,stddev) ------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; if verbose eq 1 then begin print, '' print, '--------------------------------------------------------------------------------' print, ' III. Build Taxonomic Classification' endif ;--III.1-- Loop Over Classes ------------------------------------------------------------; 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 kC=0, nbClass-1 do begin ;--III.2--Select Only Asteroids from Current Class -----------------------------------; locClass = where( strCmp(ref.class, class[kC].label), nbLoc) print, class[kC].label+'-types ('+strtrim(string(nbLoc),2)+')' ;--III.3--Compute Correlations and Covariances ---------------------------------------; locRef = ref[locClass].col corr = correlate(locRef) cov = correlate(locRef, /cov ) ;--III.3--Determine Class Orientation and Deviation ----------------------------------; for kC1=0, nbComb-1 do begin ;--III.3.1--Sample Variance & Normalization of Covariance varComb1 = variance( ref[locClass].col[kC1] ) cov[kC1,*] /= varComb1 ;--III.3.2--Loop over Second Dimension for kC2=kC1+1, nbComb-1 do begin ;--III.3.3--Keep Correlation and Covariance covArr[kC,kC1,kC2].correl = corr[kC1,kC2] covArr[kC,kC1,kC2].covar = cov[kC1,kC2] ;--III.3.4--Compute Cloud Center covArr[kC,kC1,kC2].meanX = mean( ref[locClass].col[kC1] ) covArr[kC,kC1,kC2].meanY = mean( ref[locClass].col[kC2] ) covArr[kC,kC1,kC2].stdX = stddev( ref[locClass].col[kC1] ) covArr[kC,kC1,kC2].stdY = stddev( ref[locClass].col[kC2] ) ;--III.3.5--Compute Cloud Orientation covArr[kC,kC1,kC2].angle = atan( covArr[kC,kC1,kC2].covar ) covArr[kC,kC1,kC2].cos = cos( covArr[kC,kC1,kC2].angle ) covArr[kC,kC1,kC2].sin = sin( covArr[kC,kC1,kC2].angle ) ;--III.3.6--Compute Cloud Dispersion newX = ref[locClass].col[kC1] - covArr[kC,kC1,kC2].meanX newY = ref[locClass].col[kC2] - covArr[kC,kC1,kC2].meanY rotX = newX*covArr[kC,kC1,kC2].cos + newY*covArr[kC,kC1,kC2].sin rotY =-newX*covArr[kC,kC1,kC2].sin + newY*covArr[kC,kC1,kC2].cos covArr[kC,kC1,kC2].stdA = stddev( rotX ) covArr[kC,kC1,kC2].stdB = stddev( rotY ) endfor endfor endfor ;--III.4-- Create Validation Figures ----------------------------------------------------; ;--III.4.1-- Single File with all plots if plotReferenceMulti eq 1 then $ movisPlotMulti, dirTaxo+'Reference.eps', class, covArr, ref, range=ccRang, label=combLabel ;--III.4.2-- Multiple Files for each plot if plotReferenceSingle eq 1 then $ for kC1=0, nbComb-1 do for kC2=nbComb-1, kC1+1, -1 do begin singleName = dirTaxo+'Reference-'+combLabel[kC1]+'_vs_'+combLabel[kC2]+'.eps' singleCovar = covArr[*,kC1,kC2] singleRange = ccRang[*,[kC1,kC2]] singleLabel=combLabel[[kC1,kC2]] emptySing = {num:0, name:'', class:'', N:0, col:fltarr(2)} singleRef = replicate(emptySing,nbRef) singleRef.class=ref.class singleRef.col = ref.col[[kC1,kC2]] movisPlotSingle, singleName, class, singleCovar, singleRef, range=singleRange, label=singleLabel endfor ;--III.5-- Weight for Color-Class Distances ---------------------------------------------; weightArr=fltarr(nbComb,nbComb)+1. ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--- TAG --- IV --- Read Data Set to be Classified ------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; if verbose eq 1 then begin print, '' print, '--------------------------------------------------------------------------------' print, ' IV. Read Observed Data Set' endif ;--IV.1-- Read MOVIS-Color Catalogue ----------------------------------------------------; if classTaxo eq 0 then begin ;--IV.1.1-- Read MOVIS Catalogue movis = vistaReadMOVIS( dirDATA+'MOVIS-C.csv' ) nbCatalog = n_elements(movis) ;--IV.1.2-- Apply Selection on Color Uncertainty selectMOVIS = where( movis.YmJ.unc le colorCut.YJ and $ movis.YmH.unc le colorCut.YH and $ movis.YmK.unc le colorCut.YK and $ movis.JmH.unc le colorCut.JH and $ movis.JmK.unc le colorCut.JK and $ movis.HmK.unc le colorCut.HK, nbSSO) movis = movis[selectMOVIS] ;--IV.1.3-- Create the Data Structure ssoCol = [[movis.YmJ.val],[movis.YmH.val],[movis.YmK.val],[movis.JmH.val],[movis.JmK.val],[movis.HmK.val]] sso = replicate(emptySSO,nbSSO) sso.num = movis.num sso.name = movis.name sso.class= replicate('U',nbSSO) for k=0,nbSSO-1 do sso[k].col = ssoCol[k,*] ;--IV.2-- Use Reference Sample as Fake Observing Sample ---------------------------------; endif else begin ;--IV.2.1-- Read Bus-DeMeo Name readcol, dirDATA+'bdm.csv', bdmNum, bdmName, bdmC1, bdmC2, bdmC3, bdmDate, $ format='(L,A,A,A,A,A)', delimiter=',', /Silent ;--IV.2.2-- Select ALL Spectra nbCatalog = n_elements(bdmNum) nbSSO = nbCatalog ;--IV.2.3-- Create the Data Structure sso = replicate(emptySSO,nbSSO) sso.num = bdmNum sso.name = bdmName sso.class= strTrim(bdmC3,2) sso.col = ref.col ;--TBD single class in ref if keyword_set(bdmChoice) then begin ;--TBD single class in ref select = where( strCmp( bdmChoice, strMid(strTrim(bdmC3,2),0,1) ), nbSSO ) ;--TBD single class in ref ssoNum = ssoNum[select] ;--TBD single class in ref ssoName = ssoName[select] ;--TBD single class in ref ssoKnown = ssoKnown[select] ;--TBD single class in ref ssoCol = ssoCol[select,*] ;--TBD single class in ref endif ;--TBD single class in ref selectMOVIS=indgen(nbSSO) ;--TBD single class in ref print, 'Number of asteroids in Bus-DeMeo:', nbSSO endelse ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; ;--- TAG --- V --- Classify SSO using Ellipses from Reference ------------------------; ;------------------------------------------------------------------------------------------; ;------------------------------------------------------------------------------------------; if verbose eq 1 then begin print, '' print, '--------------------------------------------------------------------------------' print, ' V. Class the Observed Sample' endif for kSSO=0, nbSSO-1 do begin ;--V.1-- Metric Arrays ----------------------------------------------------------------; distArr=fltarr(nbClass) ;- Distance in color-color dimension chi2Arr=fltarr(nbClass) ;- Chi2 in color-color probArr=fltarr(nbClass) ;- Probability associated with chi2 UsedFilter = intarr(4) ;- List of used filters ;-mp ;-mp ;-mp ;--IV.1--Loop Over SSO in MOVIS --------------------------------------------------------- ;-mp ssoClass = replicate('U',nbSSO) ;-mp ssoComb = intarr(nbSSO) ;-mp;struc = fl ;-mp; struc = astelements(433) ;-mp ;-mp fname='~/W/Vista/datos/resumen.dat' ; to write down results from line 712 in file resumen.dat ;+++++++++++++++++++++ ;-mp OPENW,1,fname ;+++++++++++++++++++++ ;-mp ;-mp printf, 1, 'Number', 'Name', 'Class', 'id','A','B','C','D','K','L','Q','S','T','X','V', $ ;-mp 'Y-J','Y-H','Y-Ks','J-H ','J-Ks','H-Ks','a','e','i',$ ;-mp format='(A-6,2x,A-20,6x,A-10,4x,11(A-2,6x),6(A-3,7x),3(A-1,10x))' ;-mp ;-mp; struc = astelements(kSSO) ;-mp ;--V.2-- Loop Over Taxonomic Classes --------------------------------------------------; for kC=0, nbClass-1 do begin chi2 = 0 sso[kSSO].N=0 ;--V.3--Loop Over Valid Combination of Filters --------------------------------------; for kC1=0, nbComb-1 do begin if sso[kSSO].col[kC1] ne -99.99 then $ ; for kC2=round(nbComb)-1, kC1+1, -1 do begin for kC2=kC1+1, nbComb-1 do begin if sso[kSSO].col[kC2] ne -99.99 then begin ;--V.4-- Compute Chi2 to Current Class --------------------------------------; dX = sso[kSSO].col[kC1]-covArr[kC,kC1,kC2].meanX dY = sso[kSSO].col[kC2]-covArr[kC,kC1,kC2].meanY rotX = dX*covArr[kC,kC1,kC2].cos + dY*covArr[kC,kC1,kC2].sin rotY =-dX*covArr[kC,kC1,kC2].sin + dY*covArr[kC,kC1,kC2].cos chi2 += ( (rotX/covArr[kC,kC1,kC2].stdA)^2. + $ (rotY/covArr[kC,kC1,kC2].stdB)^2. ) * weightArr[kC1,kC2] sso[kSSO].N++ ;--V.5-- Keep Track of Used Filters -----------------------------------------; ;--V.5.1-- First Axis case combLabel[kC1] 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 ;--V.5.2-- Second Axis case combLabel[kC2] 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 ;--V.5.3-- Once all four filters have been used -> exit if total(UsedFilter) eq 4 then begin kC1 = nbComb kC2 = nbComb endif endif endfor endfor ;--End of V.3--Loop Over Valid Combination of Filters ;--V.6--Degrees of Freedom, Distance, Probability -----------------------------------; DoF = TOTAL(UsedFilter)-1 chi2Arr[kC] = chi2 distArr[kC] = sqrt(chi2) probArr[kC]= 1 - chisqr_pdf(chi2, DoF) endfor ;--End of V.2-- Loop Over Taxonomic Classes ;--V.7--Assign Class: Multi-Colors --------------------------------------------------; if sso[kSSO].N ge 1 then begin ;--V.7.1-- Select Probability over Threshold selProba = where( probArr ge minProba, nbSelProba ) if nbSelProba ge 1 then begin ;--V.7.1/A-- Sort Classes in Probability preOrder = sort( chi2Arr[selProba] ) order = selProba[preOrder] ;--V.7.1/B-- Store Probability & Assign Most-Likely Class ;-TBD - Room for improvement sso[kSSO].proba = probArr sso[kSSO].class = class[order[0]].label ;--V.7.2-- Objects with Low Probability Only endif else begin ;-TBD Myriam endelse ;--V.8--Assign Class: Single-Color --------------------------------------------------; endif else begin ;--V.8.1-- What Color Was It? pair = where( sso[kSSO].col[*] ne -99.99 ) color = sso[kSSO].col[pair] ;--V.8.2-- Strong cut-off for specific classes case combLabel[pair] of 'Y-H' : if color lt 0.45 then sso[kSSO].class = 'B' else $ if color gt 1.15 then sso[kSSO].class = 'A' 'Y-J' : if color lt 0.25 then sso[kSSO].class = 'B' else $ if color gt 0.60 then sso[kSSO].class = 'V' 'Y-Ks': if color lt 0.50 then sso[kSSO].class = 'B' else $ if color gt 1.30 then sso[kSSO].class = 'A' 'J-H' : if color lt 0.20 then sso[kSSO].class = 'V' else $ if color gt 0.65 then sso[kSSO].class = 'A' 'J-Ks': if color lt 0.17 then sso[kSSO].class = 'V' else $ if color gt 0.80 then sso[kSSO].class = 'A' 'H-Ks': if color gt -0.03 then sso[kSSO].class = 'V' else $ if color gt 0.20 then sso[kSSO].class = 'D' else: stop endcase endelse endfor ;--End of V-- Loop over SSOs for kC=0, nbClass-1 do begin pC = where( strCmp(sso.class, class[kC].label) ) sso[pC].id = kC endfor ;--V.9-- Find Orbital Elements ----------------------------------------------------------; if findAEI eq 1 then begin wn = where( sso.num ne 0, nWN, comp=wo, nComp=nWO ) ;--V.9.1-- Numbered SSO for k=0, nWN-1 do begin orb = astElements(sso[wn[k]].num) sso[wn[k]].a = orb.a sso[wn[k]].e = orb.e sso[wn[k]].i = orb.i sso[wn[k]].H = orb.H endfor ;--V.9.1-- Numbered SSO for k=0, nWO-1 do begin orb = astElements(sso[wo[k]].name) sso[wo[k]].a = orb.a sso[wo[k]].e = orb.e sso[wo[k]].i = orb.i sso[wo[k]].H = orb.H endfor endif ;--V.10-- Store results on Disk ---------------------------------------------------------; ;--V.10.1--Open File & Header openW, idEx, dirDATA+'movis-classed.csv', /Get_Lun header='Number, Name, Class, classId, N, ' for kC=0, nbClass-1 do header+= class[kC].label+', ' for kC=0, nbComb-1 do header+= combLabel[kC]+', ' header += 'a, e, i, H' printf, idEx, header for kSSO=0, nbSSO-1 do begin ;--V.10.2-- Basic Information line = string(sso[kSSO].num,format='(I6)')+', '+$ string(sso[kSSO].name,format='(A-20)')+', '+$ string(sso[kSSO].class,format='(A1)')+', '+$ string(sso[kSSO].id,format='(I2)')+', '+$ string(sso[kSSO].N,format='(I2)')+', ' ;--V.10.3-- Probability for Each Class for kC=0, nbClass-1 do line+= string(sso[kSSO].proba[kC],format='(F6.4)')+', ' ;--V.10.4-- Colors from MOVIS-C for kC=0, nbComb-1 do line+= string(sso[kSSO].col[kC],format='(F8.4)')+', ' ;--V.10.5-- Orbital Elements line+= string(sso[kSSO].a, format='(D10.5)')+', '+$ string(sso[kSSO].e, format='(F9.6)')+', '+$ string(sso[kSSO].i, format='(F9.3)')+', '+$ string(sso[kSSO].H, format='(F6.2)') ;--V.10.6-- Push line to File printf, idEx, line endfor free_lun, idEX stop ; forprint, sso.num, sso.name, sso.class, $ ; sso.proba[0], sso.proba[1], sso.proba[2], sso.proba[3], sso.proba[4], $ ; sso.proba[5], sso.proba[6], sso.proba[7], sso.proba[8], sso.proba[9], sso.proba[10], $ ; 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], sso[kSSO].col[0], sso[kSSO].col[1], sso[kSSO].col[2],$ ; sso[kSSO].col[3], sso[kSSO].col[4], sso[kSSO].col[5], struc.a, struc.e,struc.i,$ ;; sso[kSSO].col[3], sso[kSSO].col[4], sso[kSSO].col[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,", "))' ; CLOSE,1 ; ============================================================================== a ver ;stop cghistoplot, ssoComb, binsize=1, histdata=histComb, loca=xComb print, 'Number of asteroids in the catalog:', nbCatalog print, 'Number after uncertainty selection:', nbSSO 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 minNumOfColor colors ------------------------------------ sel = where( ssoComb ge minNumOfColor, nbSel ) print, 'Number with more than '+strtrim(string(minNumOfColor),2)+' colors:', nbSel if checkTaxo eq 0 then begin index = where( ssoComb ge minNumOfColor 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 kC1=0, nbComb-1 do begin for kC2=kC1+1, nbComb-1 do begin ;--V.2.2--Prepare Plot OutLay cgPlot, ssoCol[index,kC1], ssoCol[index,kC2], /NoData, xTitle=combLabel[kC1], yTitle=combLabel[kC2] ;--V.2.3--Plot Each Object within the Class for kC=0, nbClass-1 do begin cur=where( strcmp(class[kC].label,ssoClass) ) cgPlot, /OverPlot, ssoCol[cur,kC1], ssoCol[cur,kC2], $ psym=class[kC].sym, color=class[kC].color, symSize=class[kC].size endfor ;--V.2.4--Plot Cross + Ellipse for each Class for kC=0, nbClass-1 do begin cgPlot, covArr[kC,kC1,kC2].meanX+covArr[kC,kC1,kC2].stdA*cos(theta), $ covArr[kC,kC1,kC2].meanY+covArr[kC,kC1,kC2].stdB*sin(theta), color=class[kC].color, /over cgErrPlot, covArr[kC,kC1,kC2].meanX, color=class[kC].color, $ covArr[kC,kC1,kC2].meanY-covArr[kC,kC1,kC2].stdB, $ covArr[kC,kC1,kC2].meanY+covArr[kC,kC1,kC2].stdB cgErrPlot, covArr[kC,kC1,kC2].meanY, color=class[kC].color, $ covArr[kC,kC1,kC2].meanX-covArr[kC,kC1,kC2].stdA, $ covArr[kC,kC1,kC2].meanX+covArr[kC,kC1,kC2].stdA, /horizontal cgText, covArr[kC,kC1,kC2].meanX+covArr[kC,kC1,kC2].stdA*0.3, $ covArr[kC,kC1,kC2].meanY+covArr[kC,kC1,kC2].stdB*0.3, class[kC].label, color=class[kC].color 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(minNumOfColor),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=dirFIG+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 kC1=0, nbComb-1 do begin for kC2=nbComb-1, kC1+1, -1 do begin ;--V.3.3--Prepare Plot OutLay ;--V.3.3/A--XY Range lonX = [ covArr[*,kC1,kC2].meanX+covArr[*,kC1,kC2].stdA, $ covArr[*,kC1,kC2].meanX-covArr[*,kC1,kC2].stdA, ssoCol[index,kC1]] lonY = [ covArr[*,kC1,kC2].meanY+covArr[*,kC1,kC2].stdB, $ covArr[*,kC1,kC2].meanY-covArr[*,kC1,kC2].stdB, ssoCol[index,kC2]] xRang = ccRang[*,kC1] yRang = ccRang[*,kC2] ;--V.3.3/B--Graph position pos= [xBot + kC1*xLen, $ yBot + (kC2-1)*xLen, $ xBot + (kC1+1)*xLen, $ yBot + (kC2 )*xLen ] ;--V.3.3/C--XY Axes On the First column if kC1 eq 0 then begin ;--Last line if kC2 eq kC1+1 then begin cgPlot, lonX, lonY, /NoData, /NoErase, /Normal, position=pos, $ xStyle=1, xRange=xRang, xTitle=combLabel[kC1], $ yStyle=1, yRange=yRang, yTitle=combLabel[kC2] ;--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=combLabel[kC2] endelse ;--V.3.3/D--XY Axes On the other columns endif else begin ;--Last line if kC2 eq kC1+1 then begin cgPlot, lonX, lonY, /NoData, /NoErase, /Normal, position=pos, $ xStyle=1, xRange=xRang, xTitle=combLabel[kC1], $ 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,kC1], ssoCol[cur,kC2], $ ; 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 kC=0, nbClass-1 do begin cur=where( strcmp(class[kC].label,ssoClass) and ssoComb ge minNumOfColor, nbCur ) ; print, combLabel[kC1], combLabel[kC2], class[kC].label, nbCur, format='(A-5,A-5,A-5,I5)' cgPlot, /OverPlot, ssoCol[cur,kC1], ssoCol[cur,kC2], $ psym=class[kC].sym, color=class[kC].color, symSize=class[kC].size endfor ;--V.3.5--Plot Cross + Ellipse for each Class for kC=0, nbClass-1 do begin ;--V.3.5/A--Long-Axis Line xA = [-1,1]*covArr[kC,kC1,kC2].stdA yA = [0,0] rotXA = xA*covArr[kC,kC1,kC2].cos - yA*covArr[kC,kC1,kC2].sin + covArr[kC,kC1,kC2].meanX rotYA = xA*covArr[kC,kC1,kC2].sin + yA*covArr[kC,kC1,kC2].cos + covArr[kC,kC1,kC2].meanY cgPlot, rotXA, rotYA, color=class[kC].color, /OverPlot ;--V.3.5/B--Short-Axis Line xB = [0,0] yB = [-1,1]*covArr[kC,kC1,kC2].stdB rotXB = xB*covArr[kC,kC1,kC2].cos - yB*covArr[kC,kC1,kC2].sin + covArr[kC,kC1,kC2].meanX rotYB = xB*covArr[kC,kC1,kC2].sin + yB*covArr[kC,kC1,kC2].cos + covArr[kC,kC1,kC2].meanY cgPlot, rotXB, rotYB, color=class[kC].color, /OverPlot ;--V.3.5/C--Ellipse x = covArr[kC,kC1,kC2].stdA * cos(theta) y = covArr[kC,kC1,kC2].stdB * sin(theta) rotX = x*covArr[kC,kC1,kC2].cos - y*covArr[kC,kC1,kC2].sin + covArr[kC,kC1,kC2].meanX rotY = x*covArr[kC,kC1,kC2].sin + y*covArr[kC,kC1,kC2].cos + covArr[kC,kC1,kC2].meanY cgPlot, rotX, rotY, color=class[kC].color, /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 kC=0, nbClass-1 do begin cur=where( strcmp(class[kC].label,ssoClass) and ssoComb ge minNumOfColor, nbCur ) ; cur=where( strcmp(ssoClass,class[kC].label), nbCur ) cgPlot, xKey, yKey+kC*yShi, psym=class[kC].sym, color=class[kC].color, symSize=class[kC].size, /OverPlot cgText, xKey+xShi, yKey+kC*yShi+yTxt, class[kC].label+'-types: ', color=class[kC].color cgText, xKey+xShi*13, yKey+kC*yShi+yTxt, strtrim(string(nbCur),2), align=1, color=class[kC].color endfor cgPS_close, /png, Delete_PS=0 ;--V.4--Validation Figures in Many Files ---------------------------------------------- ;--V.4.1-- Loop over Combinations for kC1=0, nbComb-1 do begin for kC2=nbComb-1, kC1+1, -1 do begin ;--V.4.2-- Prepare plot if checkTaxo eq 0 then begin filename = 'MOVIS-Result-'+strtrim(string(minNumOfColor),2)+'-'+strtrim(string(medCut,format='(F4.2)'),2)+$ '-'+combLabel[kC1]+'_vs_'+combLabel[kC2]+'.eps' endif else begin fileName = 'DeMeo-Valid-'+combLabel[kC1]+'_vs_'+combLabel[kC2]+'.eps' endelse cgPS_open, Filename=dirFIG+fileName, /metric, /decomposed, /encapsulated, $ xSize=15, ySize=15, language_level=2, /quiet ;--V.4.3--Prepare Plot OutLay xRang = ccRang[*,kC1] yRang = ccRang[*,kC2] cgPlot, 0,0, /NoData, charsize=1, $ xStyle=1, xRange=xRang, xTitle=combLabel[kC1], $ yStyle=1, yRange=yRang, yTitle=combLabel[kC2], $ xTickInt=ccTickInt[kC1], xMinor=ccTickMin[kC1], $ yTickInt=ccTickInt[kC2], yMinor=ccTickMin[kC2] ; cur=where( strcmp('U',ssoClass), nbCur ) ; ++++++++++++++++++++++++= cgPlot, /OverPlot, ssoCol[cur,kC1], ssoCol[cur,kC2], $ psym=cSymb[nbClass], color=cColor[nbClass], symSize=sSymb[nbClass] ;++++++++++++++++++ def symb col size ;--V.4.4--Plot Each Object within the Class for kC=0, nbClass-1 do begin cur=where( strcmp(class[kC].label,ssoClass) and ssoComb ge minNumOfColor, nbCur ) ; cur=where( strcmp(class[kC].label,ssoClass), nbCur ) ; print, combLabel[kC1], combLabel[kC2], class[kC].label, nbCur, format='(A-5,A-5,A-5,I5)' cgPlot, /OverPlot, ssoCol[cur,kC1], ssoCol[cur,kC2], $ psym=class[kC].sym, color=class[kC].color, symSize=class[kC].size endfor ;--V.4.5--Plot Cross + Ellipse for each Class for kC=0, nbClass-1 do begin ;--V.4.5/A--Long-Axis Line xA = [-1,1]*covArr[kC,kC1,kC2].stdA yA = [0,0] rotXA = xA*covArr[kC,kC1,kC2].cos - yA*covArr[kC,kC1,kC2].sin + covArr[kC,kC1,kC2].meanX rotYA = xA*covArr[kC,kC1,kC2].sin + yA*covArr[kC,kC1,kC2].cos + covArr[kC,kC1,kC2].meanY cgPlot, rotXA, rotYA, color=class[kC].color, /OverPlot ;--V.4.5/B--Short-Axis Line xB = [0,0] yB = [-1,1]*covArr[kC,kC1,kC2].stdB rotXB = xB*covArr[kC,kC1,kC2].cos - yB*covArr[kC,kC1,kC2].sin + covArr[kC,kC1,kC2].meanX rotYB = xB*covArr[kC,kC1,kC2].sin + yB*covArr[kC,kC1,kC2].cos + covArr[kC,kC1,kC2].meanY cgPlot, rotXB, rotYB, color=class[kC].color, /OverPlot ;--V.4.5/C--Ellipse x = covArr[kC,kC1,kC2].stdA * cos(theta) y = covArr[kC,kC1,kC2].stdB * sin(theta) rotX = x*covArr[kC,kC1,kC2].cos - y*covArr[kC,kC1,kC2].sin + covArr[kC,kC1,kC2].meanX rotY = x*covArr[kC,kC1,kC2].sin + y*covArr[kC,kC1,kC2].cos + covArr[kC,kC1,kC2].meanY cgPlot, rotX, rotY, color=class[kC].color, /OverPlot cgText, covArr[kC,kC1,kC2].meanX + 0.3*covArr[kC,kC1,kC2].stdX, $ covArr[kC,kC1,kC2].meanY + 0.3*covArr[kC,kC1,kC2].stdY, class[kC].label, color=class[kC].color 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 kC=0, nbClass-1 do begin cur = where( strcmp( Class, class[kC].label), 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( class[kC].label, 'T' ) then begin ord=sort(numPerClass) endif xmN[kC,0:nbUniq-1] = numPerClass[ord] xmF[kC,0:nbUniq-1] = frqPerClass[ord] xmL[kC,0:nbUniq-1] = uniqClass[ord] print, '' print, class[kC].label, ' -- ', 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=dirFIG+'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(class.Label), 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 kC=nbClass-1,0,-1 do begin for kC=0, nbClass-1 do begin p=where( strcmp(class.label, xmL[kC,kE]) ) if p[0] ne -1 then begin colorVec[kC] = cColor[p[0]] if comArr[kC,p[0]] lt 2 then xmV[kC,kE]=1 statValid[*, comArr[kC,p[0]]] += xmF[*,kC] if xmV[kC,kE] eq 1 and xmF[kC,kE] gt 0.5 then begin base = last last[kC] += xmF[kC,kE] xx = base[kC] + [0,xmF[kC,kE],xmF[kC,kE],0,0] yy = [0,0,1,1,0]*barSize + barOffSet + yArr[nbClass-1-kC] cgPlot, /over, xx, yy, color=colorVec[kC] cgText, mean(xx[0:1]), yy[1] + (yy[2]-yy[1])/4., align=0.5, xmL[kC,kE], color=colorVec[kC] 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