$Id$
@read_sc_1d_table
@read_sc_2d_table
@read_sc_3d_table
@get_sc_prob_1d
@get_sc_prob_2d
@get_sc_prob_3d
@get_sc_prior_1d
@get_sc_prior_2d
@get_sc_prior_3d
@get_sc_prior_map
@compute_month_day
@leapyr
;--------------------------------------------------------------------------------------
; test IDL implentation of new naive bayes cloud mask/phase
;--------------------------------------------------------------------------------------
pro test_new_sc, training_file = training_file, nskip = nskip

if (~keyword_set(nskip)) then nskip = 10

;--------------------------------------------------------------------------------------
;-- constants
;--------------------------------------------------------------------------------------
missing = -999.0
cloud_clear_prob_thresh = 0.5
prob_min_thresh = 0.001
r_max = 1000.0
cod_clear_thresh = 0.100000
cod_ice_water_thresh = 0.100000
lat_min = -90.0
lat_max = 90.0
etrop_cloud_thresh = 0.9
count_min = 10.0
nsfc = 7

solzen_min = 0.0
solzen_max = 180.0
glintzen_min = 0.0
glintzen_max = 180.0
zsfc_min = -100.0
zsfc_max = 20000.0
snow_class_min = 1
snow_class_max = 3

;--------------------------------------------------------------------------------------
;--- read prior
;--------------------------------------------------------------------------------------
sc_prior_map_name = 'sc_prob_calipso_prior.nc'
get_sc_prior_map, sc_prior_map_name

;--------------------------------------------------------------------------------------
;--- read in observations
;--------------------------------------------------------------------------------------
restore, training_file

nx = n_elements(lat)

logcod065 = cod065
idx = where(cod065 gt 0.0,cc) & if (cc gt 0) then logcod065[idx] = alog10(cod065[idx])

logcod138 = cod138
idx = where(cod138 gt 0.0,cc) & if (cc gt 0) then logcod138[idx] = alog10(cod138[idx])

;-------------------------------------------------------------------------------------
;--- filter lidar
;-------------------------------------------------------------------------------------
;--- very thin goes to clear
idx = where(cod_lidar ne missing and cod_lidar lt cod_clear_thresh,cc) & if (cc gt 0) then phase_lidar[idx] = 0
;--- very thin on top of water goes to water
idx = where(cod_lidar gt cod_clear_thresh and codi_lidar ne missing and codi_lidar lt cod_ice_water_thresh,cc) & if (cc gt 0) then phase_lidar[idx] = 2
;--- partial cloudiness is set to bad
idx = where(frac_lidar ne 0.0 and frac_lidar ne 1.0,cc) & if (cc gt 0) then phase_lidar[idx] = -1

;--------------------------------------------------------------------------------------
;-- define sensor
;--------------------------------------------------------------------------------------
sensor = 'modis'

;--------------------------------------------------------------------------------------
;-- define table
;--------------------------------------------------------------------------------------
@sc_table_library.inc

;--------------------------------------------------------------------------------------
;--- read tables
;-------------------------------------------------------------------------------------
table_read_flag[*] = 0
for itable = 0, ntables - 1 do begin

    temp_table_name = sensor + "_" + table_name[itable] + "_sc_prob.nc"

    if (table_rank[itable] eq 1) then begin
       read_sc_1d_table, temp_table_name, lut_temp, attrs_temp, read_flag
       *(lut[itable]) = lut_temp
       *(attrs[itable]) = attrs_temp
       table_read_flag[itable] = read_flag
    endif

    if (table_rank[itable] eq 2) then begin
       read_sc_2d_table, temp_table_name, lut_temp, attrs_temp, read_flag
       *(lut[itable]) = lut_temp
       *(attrs[itable]) = attrs_temp
       table_read_flag[itable] = read_flag
    endif

    if (table_rank[itable] eq 3) then begin
       read_sc_3d_table, temp_table_name, lut_temp, attrs_temp, read_flag
       *(lut[itable]) = lut_temp
       *(attrs[itable]) = attrs_temp
       table_read_flag[itable] = read_flag
    endif

    ;--- turn off tables that were not read in
    if (table_read_flag[itable] eq 0) then table_on_flag[itable] = 0

endfor

;-------------------------------------------------------------------------------------
; check that the channels are on for each table
;-------------------------------------------------------------------------------------
for itable = 0, ntables - 1 do begin
    for i = 0, nchan_tables - 1 do begin
      if(chan_on_table[itable,i] gt 0) then begin
         if(chan_on_flag[chan_on_table[itable,i]] eq 0) then table_on_flag[itable] = 0
      endif
    endfor
endfor

;-------------------------------------------------------------------------------------
; define output
;-------------------------------------------------------------------------------------
post_prob_sc = make_array(nx,/FLOAT,VALUE=missing)
post_prob_nosc = make_array(nx,/FLOAT,VALUE=missing)
sc_mask_new = make_array(nx,/INT,VALUE=-1)
sc_mask_true = make_array(nx,/INT,VALUE=-1)


;------------------------------------------------------------------------------------
; define truth
;------------------------------------------------------------------------------------
idx = where(phase_lidar eq 0 or phase_lidar eq 1 or phase_lidar eq 3 or tc_lidar gt 273.15, cc)
if (cc gt 0) then sc_mask_true[idx] = 0
idx = where(phase_lidar eq 2 and tc_lidar gt 0.0 and tc_lidar le 273.15, cc)
if (cc gt 0) then sc_mask_true[idx] = 1


;-----------------------------------------------------------
; determine first valid table index
;-----------------------------------------------------------
first_valid_table = -1 
itable = 0
while (first_valid_table lt 0 and itable lt ntables) do begin
  if (table_read_flag[itable] eq 1) then first_valid_table = itable
  itable++
endwhile

if (first_valid_table lt 0) then begin
   print, 'ERROR, no valid tables found'
endif

;-----------------------------------------------------------
; Loop over Pixels and generate probabilities
;-----------------------------------------------------------
for i = 0, nx-1, nskip do begin

    isfc = bayes_mask_type[i] - 1 ; idl is zero based

    ;--- prior
    compute_month_day, fix(doy[i]), leapyr(fix(year[i])), month, day
    imonth = month - 1
    idiurnal = 0
    get_sc_prior_values, lon[i], lat[i], imonth, idiurnal, prior_prob_sc

    ;--- take prior for this from any of the luts
;   f = *(lut[first_valid_table])
;   prior_prob_cloud = f[isfc].cloud_fraction
;   prior_prob_clear =  1.0 - prior_prob_cloud
;   prior_prob_ice = f[isfc].ice_fraction 
;   prior_prob_water = f[isfc].water_fraction 

;   ;-- base prior on etrop
;   x = etrop[i]
;   get_prior, x, isfc, etrop_btd8511_lut, etrop_btd8511_attrs, prior_prob_clear, prior_prob_ice, prior_prob_water

    ;--- force water and ice priors to be equal (combats over-icing)
;   prior_prob_water = 0.5*(1.0-prior_prob_clear)
;   prior_prob_ice = prior_prob_water

    ;---- reset table probabilities to missing
    sc_prob[*] = missing
    obs_prob[*] = missing

    ;---- loop through each table
    for itable = 0, ntables - 1 do begin

        if (table_on_flag[itable] eq 0) then goto, skip_table

        x = missing
        y = missing
        z = missing

        case table_name[itable] OF

           ;--- 1d
           'all_topa': begin
              x = topa[i]
           end

           ;--- 2d
           'day_topa_ref160': begin
              x = topa[i]
              y = ref160[i]
           end
           'day_topa_ref375': begin
              x = topa[i]
              y = ref375[i]
           end
           'night_topa_ref375': begin
              x = topa[i]
              y = ref375[i]
           end
           'all_topa_btd8511': begin
              x = topa[i]
              y = t85[i] - t11[i]
           end
           'all_topa_btd1112': begin
              x = topa[i]
              y = t11[i] - t12[i]
           end

           ;--- 3d

        endcase

        if (table_rank[itable] eq 1) then begin
            get_sc_prob_1d,  x, solzen[i], glintzen[i], zsfc[i], snow_class[i], isfc,$
                          *(lut[itable]), *(attrs[itable]), prob_min_thresh, r_max, missing, $
                          z1, z2
        endif

        if (table_rank[itable] eq 2) then begin
            get_sc_prob_2d,  x, y, solzen[i], glintzen[i], zsfc[i], snow_class[i], isfc,$
                          *(lut[itable]), *(attrs[itable]), prob_min_thresh, r_max, missing, $
                          z1, z2
        endif

        if (table_rank[itable] eq 3) then begin
            get_sc_prob_3d,  x, y, z, solzen[i], glintzen[i], zsfc[i], snow_class[i], isfc,$
                          *(lut[itable]), *(attrs[itable]), prob_min_thresh, r_max, missing, $
                          z1, z2
        endif

        sc_prob[itable] = z1
        obs_prob[itable] = z2

       skip_table:

    endfor

    ;--------------------------------------------------------------------------------------------
    ;  combine probabilities from each tables into a final value using standard naive bayes     
    ;--------------------------------------------------------------------------------------------
    post_prob_sc[i] = prior_prob_sc
    post_prob_nosc[i] = 1.0 - prior_prob_sc
    for itable = 0, ntables - 1 do begin
       if (sc_prob[itable] ne missing) then post_prob_sc[i] = post_prob_sc[i] * sc_prob[itable]
       if (sc_prob[itable] ne missing) then post_prob_nosc[i] = post_prob_nosc[i] * (1.0-sc_prob[itable])
    endfor

    ;--- normalize these probabilties so that they sum to 1
    denominator = post_prob_sc[i] + post_prob_nosc[i]
    post_prob_sc[i] =  post_prob_sc[i] / denominator
    post_prob_nosc[i] =  post_prob_nosc[i] / denominator

    ;--- each test
;   for itable = 0, ntables - 1 do begin
;      post_prob_sc_test[itable,i] = prior_prob_sc   ; default
;      if (sc_prob[itable] ne missing) then begin
;           post_prob_sc_test[itable,i] = prior_prob_sc*sc_prob[itable] / $
;                      (prior_prob_sc*sc_prob[itable] + (1.0-prior_prob_sc)*(1.0-sc_prob[itable]))
;      endif
;   endfor

    ;-- construct a 2-level supercooled mask
    if (post_prob_sc[i] ge 0.0) then sc_mask_new[i] = 0
    if (post_prob_sc[i] ge 0.5) then sc_mask_new[i] = 1

endfor

;-----------------------------------------------------------------------
; destroy pointers
;-----------------------------------------------------------------------
lut = !null
attrs = !null

;----------------------------------------------------------------------
; performance
; note, could just do supercooled phase as aux
;---------------------------------------------------------------------
for isfc = 0, nsfc - 1 do begin

  idx_sfc = where(bayes_mask_type eq isfc + 1, cc)

  idx_valid = where(sc_mask_true[idx_sfc] ge 0 and sc_mask_new[idx_sfc] ge 0, count_valid)
  idx_agree_sc = where(sc_mask_true[idx_sfc[idx_valid]] eq 1 and sc_mask_new[idx_sfc[idx_valid]] eq 1, count_agree_sc)
  idx_false_sc = where(sc_mask_true[idx_sfc[idx_valid]] eq 0 and sc_mask_new[idx_sfc[idx_valid]] eq 1, count_false_sc)
  idx_missed_sc = where(sc_mask_true[idx_sfc[idx_valid]] eq 1 and sc_mask_new[idx_sfc[idx_valid]] eq 0, count_missed_sc)
  idx_agree_nosc = where(sc_mask_true[idx_sfc[idx_valid]] eq 0 and sc_mask_new[idx_sfc[idx_valid]] eq 0, count_agree_nosc)
  idx_sc_new = where(sc_mask_new[idx_sfc[idx_valid]] ge 1, count_sc_new)
  idx_sc_true = where(sc_mask_true[idx_sfc[idx_valid]] ge 1, count_sc_true)

  pc_sc_new =  (float(count_agree_nosc) + float(count_agree_sc)) / float(count_valid)
  false_sc_new =  float(count_false_sc) / float(count_valid)
  missed_sc_new =  float(count_missed_sc) / float(count_valid)
  frac_sc_new =  float(count_sc_new) / float(count_valid)
  frac_sc_true =  float(count_sc_true) / float(count_valid)
  print, "new = ", isfc, frac_sc_true, frac_sc_new, pc_sc_new, false_sc_new, missed_sc_new

  print, "--------------"

endfor

end