import numpy as np
from probsev_functions import *
def probsev_demo():
  #==sat probs set to missing by default
  prob_emiss_sev = 1.0
  prob_emiss_weak = 1.0
  prob_icecf_sev = 1.0
  prob_icecf_weak = 1.0
  
  mucape = 2506 #J/kg
  ebs = 20 #m/s
  max_mesh = 0.68 #inches
  flash_rate = 55 #flashes/min
  maxrc_emiss = 2.1 #This is the max in the lifetime of the storm (up to it's current time)!
  emiss_delta = 7 #min
  maxrc_icecf = -999 #missing. This is the max in the lifetime of the storm (up to it's current time)!
  icecf_delta = -999 #missing
  
  probs_info, bins_info, climo_probs = lookup_tables_edit('/home/jcintineo/probsevere/src/static/')
  p = probs_info
  b = bins_info

  #apriori factor
  severe_climo = compute_apriori(mucape,ebs,b['mucape_vals'],b['ebshear_vals'],climo_probs)
  weak_climo = 1 - severe_climo
  #satellite growth rates factors
  if(maxrc_emiss > 0):
    if(emiss_delta > 10):
      prob_emiss_sev,prob_emiss_weak = get_probs(maxrc_emiss,b['minval_emiss_RO_wk'],b['minval_emiss_RO_st'],\
                       b['binsize_emiss_RO_wk'],b['binsize_emiss_RO_st'],p['rc_emiss_RO_weak_probs'],p['rc_emiss_RO_strong_probs'])
    else:
      prob_emiss_sev,prob_emiss_weak = get_probs(maxrc_emiss,b['minval_emiss_RSO_wk'],b['minval_emiss_RSO_st'],\
                       b['binsize_emiss_RSO_wk'],b['binsize_emiss_RSO_st'],p['rc_emiss_RSO_weak_probs'],p['rc_emiss_RSO_strong_probs'])
  if(maxrc_icecf > 0):
    if(icecf_delta > 10):
      prob_icecf_sev,prob_icecf_weak = get_probs(maxrc_icecf,b['minval_icecf_RO_wk'],b['minval_icecf_RO_st'],\
                       b['binsize_icecf_RO_wk'],b['binsize_icecf_RO_st'],p['rc_icecf_RO_weak_probs'],p['rc_icecf_RO_strong_probs'])
    else:
      prob_icecf_sev,prob_icecf_weak = get_probs(maxrc_icecf,b['minval_icecf_RSO_wk'],b['minval_icecf_RSO_st'],\
                       b['binsize_icecf_RSO_wk'],b['binsize_icecf_RSO_st'],p['rc_icecf_RSO_weak_probs'],p['rc_icecf_RSO_strong_probs'])
  #mesh factors
  prob_mesh_sev,prob_mesh_weak = get_probs(max_mesh,b['minval_mesh_wk'],b['minval_mesh_st'],b['binsize_mesh_wk'],b['binsize_mesh_st'],\
                                   p['max_mesh_weak_probs'],p['max_mesh_strong_probs'])
  #lighting/EBS factors
  ns_fr_ebs_prob = get_2d_prob(p['ns_fr_ebs'],b['fr_2d_vals'],b['ebs_2d_vals'],flash_rate,ebs)
  sev_fr_ebs_prob = get_2d_prob(p['sev_fr_ebs'],b['fr_2d_vals'],b['ebs_2d_vals'],flash_rate,ebs)
  if(sev_fr_ebs_prob/ns_fr_ebs_prob > 20):
    sev_fr_ebs_prob = 0.1
    ns_fr_ebs_prob = 0.1/20.0 #==this is so the ratio of probs isn't crazy, like in our poorly sampled part of the joint distribution

  #----------------------------------------#      
  #---now compute the actual probability---#
  #----------------------------------------#
  prob1 = ((severe_climo)*(prob_emiss_sev)*(prob_icecf_sev)*(sev_fr_ebs_prob)*(prob_mesh_sev)) / \
          ((severe_climo)*(prob_emiss_sev)*(prob_icecf_sev)*(sev_fr_ebs_prob)*(prob_mesh_sev) + (weak_climo)*(prob_emiss_weak)*(prob_icecf_weak)*(ns_fr_ebs_prob)*(prob_mesh_weak))
  prob2 = ((severe_climo)*(prob_emiss_sev)*(sev_fr_ebs_prob)*(prob_mesh_sev)) / ((severe_climo)*(prob_emiss_sev)*(sev_fr_ebs_prob)*(prob_mesh_sev) + (weak_climo)*(prob_emiss_weak)*(ns_fr_ebs_prob)*(prob_mesh_weak))
  prob3 = ((severe_climo)*(prob_icecf_sev)*(sev_fr_ebs_prob)*(prob_mesh_sev)) / ((severe_climo)*(prob_icecf_sev)*(sev_fr_ebs_prob)*(prob_mesh_sev) + (weak_climo)*(prob_icecf_weak)*(ns_fr_ebs_prob)*(prob_mesh_weak))
  prob4 = ((severe_climo)*(sev_fr_ebs_prob)*(prob_mesh_sev)) / ((severe_climo)*(sev_fr_ebs_prob)*(prob_mesh_sev) + (weak_climo)*(ns_fr_ebs_prob)*(prob_mesh_weak))

  #We take the max so that weak growth rates in large, mature satellite objects don't hurt us.
  return max([prob1,prob2,prob3,prob4]) #==take the max prob of NWP+radar+sat+ltg and NWP+radar+ltg


if __name__ == '__main__':
  print('Probability of severe: ' + str(probsev_demo()))