function [] = machineLearnProfile(plot_importance,predictions,statistics)
%The following code extracts and calculates predictors for the dropsonde ML algorithm 

%Inputs:    campaign: name of the field campaign for machine learning
%           split: switch for randomizing files for testing (0 or 1)
%           plot_importance: switch for generating importance plot (0 or 1)
%           predictions: switch for outputting flight predictions as .txt
%           file
%           statistics: switch for outputting stats as .txt file

%Author: Julianna Christopoulos (juliechristo@g.ucla.edu)

%Last Updated: October 14, 2024

%Output:
%           Produces predictor importance plot for regression tree 
%           Writes MLH predictions to .txt files for each flight
%           Writes statistics to .txt files for each flight 
%------------------------------ Code Begin------------------------------
    % -- Define Field Campaigns to Inlcude --
    campaigns = {'DAQ_Texas','DAQ_California'};
    plot_importance = 1;
    predictions = 0;
    statistics = 0;

   % -- Gather Files and Split into Training and Testing
    numCampaigns = length(campaigns);

     MLH_filenames = {};
     MLH_pathnames = {};

    for i = 1:numCampaigns
        campaignName = campaigns{i};
        mlh_directory = sprintf('/data2/juliechristo/%s/PBLH/', campaignName);
        MLH_directories.(['direc' num2str(i)]) = mlh_directory;
        % Get list of files in the directory
        mlh_files = dir(fullfile(mlh_directory));
        % Extract filenames and store in the array
        filenames = {mlh_files.name};
        % Filter files based on extension (.ict or .h5)
        valid_extensions = {'.csv'};
        valid_files = filenames(endsWith(filenames, valid_extensions));
        % Store filtered filenames
        MLH_filenames{i} = valid_files;
        % Construct full pathnames and store in the array
        pathnames = cellfun(@(filename) fullfile(mlh_directory, filename), valid_files, 'UniformOutput', false);
        MLH_pathnames{i} = pathnames;
    end
    all_files = horzcat(MLH_filenames{:});
    [trainFiles, testFiles] = split_files(all_files, 0.8);


     % -- Make a Strucutre of All the Data Directories -- 
    names = struct();
    for i = 1:numCampaigns
        campaignName = campaigns{i};
        %-- Initialize Data Directories --
        names.(['Direc' num2str(i)]) = sprintf('/data2/juliechristo/%s/', campaignName);
    end

    %--Extract full directories for training --
    drop_train = gatherFiles(trainFiles,names,'PBLH')';

    %--Extract full directories for testing --
    drop_test = gatherFiles(testFiles,names,'PBLH')';




     %-------------------------------------Training-----------------------------
 
     %-- Get Other Variables --
     [time_tr,drop_tr] = gather_Drop(drop_train);
     indices_tr = find_index(drop_train,time_tr);
     gps_time_tr = get_gpstime(drop_train,indices_tr,'/Archive_ml/');
     [lat_tr] = gather_NavData(drop_train,indices_tr, '/gps_lat');
     [lon_tr] = gather_NavData(drop_train,indices_tr, '/gps_lon');
     lidar_tr = matchVars_MLH(drop_train,indices_tr,'/Archive/');
     [depol_tr] = matchVars_Depol(drop_train,indices_tr,'/depol_ml/');
     [p1_tr] = matchVars_Depol(drop_train,indices_tr,'/p1_ml/');
     [p2_tr] = matchVars_Depol(drop_train,indices_tr,'/p2_ml/');
     [p3_tr] = matchVars_Depol(drop_train,indices_tr,'/p3_ml/');
     [p4_tr] = matchVars_Depol(drop_train,indices_tr,'/p4_ml/');
     h_predic_tr = {p1_tr, p2_tr, p3_tr, p4_tr};     
 
     %--Compute Zenith and Solar Hour Angles--
     [zenith_tr,hangle_tr] = calcAnglen(trainFiles,gps_time_tr, lat_tr, lon_tr);
  
     %-- Initialize Variances --
     %-- Calculate for 532nm backscatter --
     [horz_variances, mlh_means] = cellfun(@(x) calcVar_MLH(x, time_tr, 10), h_predic_tr, 'UniformOutput', false);
     [hvar_p1_tr, hvar_p2_tr, hvar_p3_tr, hvar_p4_tr] = horz_variances{:};
    
     %-- Calculate for 1064nm Aerosol Depolarization Ratio --
     [hvar_depol_tr,mean_depol_tr] = calcVar_MLH(depol_tr, time_tr, 10);
      
     % -- Add on Variance Predictors to Original Height Predictors
     predic_all_tr = {p1_tr,p2_tr,p3_tr,p4_tr,hvar_p1_tr,hvar_p2_tr,hvar_p3_tr,hvar_p4_tr,hangle_tr,depol_tr};
 
     %-- Find Minimums and Maximums -- 
     [min_tr, max_tr] = cellfun(@find_range, predic_all_tr, 'UniformOutput', false);
     [min_p1,min_p2,min_p3,min_p4,min_hvar_p1,min_hvar_p2,min_hvar_p3,min_hvar_p4,min_hangle,min_depol] = min_tr{:};
 
     [max_p1,max_p2,max_p3,max_p4,max_hvar_p1,max_hvar_p2,max_hvar_p3,max_hvar_p4,max_hangle,max_depol] = max_tr{:};

     
     %Archived Mins and Maxs
     drop_tr = cellfun(@(x) x(:)', drop_tr, 'UniformOutput', false);
     [min_drop, max_drop] = find_range(drop_tr);
 
     %-- Standardize Values for Machine Learning --
     std_values_tr = cellfun(@(x, min_val, max_val) standardizeData(x, min_val, max_val), predic_all_tr, min_tr, max_tr, 'UniformOutput', false);
     [sp1_tr,sp2_tr,sp3_tr,sp4_tr,shvar_p1_tr,shvar_p2_tr,shvar_p3_tr,shvar_p4_tr,shangle_tr,sdepol_tr] = std_values_tr{:};
     [sdrop_tr] = standardizeData(drop_tr,min_drop,max_drop);
 
     %-------------------------------------Testing-----------------------------
     %-- Get Other Variables --
     [time_tst,drop_tst] = gather_Drop(drop_test);
     indices_tst = find_index(drop_test,time_tst);
     [lat_tst] = gather_NavData(drop_test,indices_tst, '/gps_lat');
     [lon_tst] = gather_NavData(drop_test,indices_tst, '/gps_lon');
     lidar_tst = matchVars_MLH(drop_test,indices_tst,'/Archive/');
     gps_time_tst = get_gpstime(drop_test,indices_tst,'/Archive_ml/');
     [depol_tst] = matchVars_Depol(drop_test,indices_tst,'/depol_ml/');
     [p1_tst] = matchVars_Depol(drop_test,indices_tst,'/p1_ml/');
     [p2_tst] = matchVars_Depol(drop_test,indices_tst,'/p2_ml/');
     [p3_tst] = matchVars_Depol(drop_test,indices_tst,'/p3_ml/');
     [p4_tst] = matchVars_Depol(drop_test,indices_tst,'/p4_ml/');
     h_predic_tst = {p1_tst, p2_tst, p3_tst, p4_tst};

     %--Compute Zenith and Solar Hour Angles--
     [zenith_tst,hangle_tst] = calcAnglen(testFiles, gps_time_tst, lat_tst, lon_tst);
 
     %-- Initialize Variances --
     %-- Calculate for 532nm backscatter --
     [horz_variances_tst,mlh_means_tst] = cellfun(@(x) calcVar_MLH(x, time_tst, 10), h_predic_tst, 'UniformOutput', false);
     [hvar_p1_tst, hvar_p2_tst, hvar_p3_tst, hvar_p4_tst] = horz_variances_tst{:};
 
     %-- Calculate for 1064nm Aerosol Depolarization Ratio --
     [hvar_depol_tst,mean_depol_tst] = calcVar_MLH(depol_tst, time_tst, 10);
 
 
     % -- Add on Variance Predictors to Original Height Predictors
     predic_all_tst = {p1_tst,p2_tst,p3_tst,p4_tst,hvar_p1_tst,hvar_p2_tst,hvar_p3_tst,hvar_p4_tst,hangle_tst,depol_tst}; 
  
     %-- Standardize Values for Machine Learning --
     std_values_tst = cellfun(@(x, min_val, max_val) standardizeData(x, min_val, max_val), predic_all_tst, min_tr, max_tr, 'UniformOutput', false);
     [sp1_tst,sp2_tst,sp3_tst,sp4_tst,shvar_p1_tst,shvar_p2_tst,shvar_p3_tst,shvar_p4_tst,shangle_tst,sdepol_tst] = std_values_tst{:};
     [sdrop_tst] = standardizeData(drop_test,min_drop,max_drop);
%  
%     % -- Put the training and testing predictros into a table and output as
%     % a .mat -- 
     %Flatten the Training Arrays 
     train_drop = cat(2, sdrop_tr{:})';
     % Initialize flattened array
     train_data = [];
     % Flatten vertically along columns
     for i = 1:numel(std_values_tr)
         train_data = [train_data; cat(2, std_values_tr{i}{:})];
     end
     %Write into Columns
     train_data = train_data';
 
       %Edit the Columns for training linear regression and neural (no nans)   
      medians_train = zeros(1, size(train_data, 2));
      train_data_clean = train_data;
      for col = 1:size(train_data_clean,2)
          columnData = train_data_clean(:,col);
          nanIndices = isnan(columnData);
          disp(size(nanIndices))
          median_value = nanmedian(columnData);
          medians_train(col) = median_value;
          columnData(nanIndices) = median_value;
          train_data_clean(:,col) = columnData;
      end
  
      %Do the same for the responses
      train_drop_clean = train_drop;
      medianDrop = nanmedian(train_drop_clean);
      nanDrop = isnan(train_drop_clean);
      train_drop_clean(nanDrop) = medianDrop;
 
 
     %------------------------Do the Machine Learning---------------------------
 
      %--Linear Regression Model--
      linear = fitrlinear(train_data_clean,train_drop_clean,'Regularization','ridge');

      %--Tree Ensembles Model--
      treeTemplate = templateTree('Surrogate','on');
      ensemble = fitrensemble(train_data,train_drop,'Method','Bag','Learners',treeTemplate);

      %--Neural Networks model--
      neural = fitrnet(train_data_clean,train_drop_clean);
 
 
     %--Predictor Importance--

    if plot_importance == 1
        % Calculate the predictor importance
        importance = predictorImportance(ensemble);
        predictor_labels = {'P1','P2','P3','P4','P5','P6','P7','P8','P9','P10'};
    
        % Plot the bar chart
        bar_handle = bar(importance);
    
        % Set the color for each bar
        bar_handle.EdgeColor = 'none';
        % Set y-axis limits
        ylim([0 5e-4]);
    
        % Adjust x-axis tick positions and labels
        xlabel('Predictor','FontWeight','bold','FontSize',12);
        ylabel('Importance','FontWeight','bold','FontSize',12);
        xticks(1:numel(importance));
        xticklabels(predictor_labels);
        xtickangle(90); % Optional: Rotate the x-axis labels if needed
        title('Land Ensemble Model','FontWeight','bold','FontSize',14);

        ax = gca;
        set(gca, 'XTickLabel', get(gca, 'XTickLabel'), 'FontWeight', 'bold'); 
        % Save the plot as an image file
        print('importance_land_mdl_bold.png', '-dpng', '-r300');
    end
 
  
      mlh_linear = cell(length(sp1_tst), 1);
      mlh_ens = cell(length(sp1_tst), 1);
      mlh_neural = cell(length(sp1_tst), 1);
      %Make Predictions for each flight
      for i = 1:length(testFiles)
          test_MLH = sp1_tst{i};
          nan_index = find(isnan(test_MLH));
          test_data = [sp1_tst{i}', sp2_tst{i}', sp3_tst{i}', sp4_tst{i}',shvar_p1_tst{i}', shvar_p2_tst{i}', shvar_p3_tst{i}', shvar_p4_tst{i}',shangle_tst{i}',sdepol_tst{i}'];
          test_data_clean = test_data;

          %Fix nans for linear and neural problem
          for column = 1:size(test_data_clean,2)
              colNans_test = isnan(test_data_clean(:, column));
              test_data_clean(colNans_test, column) = medians_train(column);
          end
 
          %Linear Prediction + recover units
          linear_dat = predict(linear,test_data_clean).* (max_drop - min_drop) + min_drop;
          linear_dat(nan_index) = nan;
          mlh_linear{i} = linear_dat';
          %Ensemble Prediction + recover units
          ens_dat = predict(ensemble,test_data).* (max_drop - min_drop) + min_drop;
          ens_dat(nan_index) = nan;
          mlh_ens{i} = ens_dat';
          %Neural Prediction + recover units
          neural_dat = predict(neural,test_data_clean).* (max_drop - min_drop) + min_drop;
          neural_dat(nan_index) = nan;
          mlh_neural{i} = neural_dat';
      end
 
 
     %---------Ouput the Predictions for each flight-------------------------- 
     if predictions == 1
         for o = 1:length(mlh_linear)
             time_col = gps_time_tst{o}';
             linear_col = mlh_linear{o}';
             ens_col = mlh_ens{o}';
             neural_col = mlh_neural{o}';
             drop_col = drop_tst{o};
             lidar_col = lidar_tst{o}';
  
             data = table(time_col,linear_col, ens_col, neural_col,drop_col,lidar_col,'VariableNames', {'Time','Linear', 'Ensemble', 'Neural','Drop','Lidar'});
             currentfile = testFiles{o};
             date = currentfile(1:8);
             filename = sprintf('%s.txt', date);
             writetable(data, filename, 'Delimiter', '\t');
         end
     end 
     %-----------------------Compute Statistics of ML Predictions----------------------
  
     %---------------------------------------Statistics for the entire campaign--------------------------------
     if statistics == 1 
         dropsonde_all = [];
         lidar_all = [];
         def_all = [];
         ensemble_all = [];
         linear_all = [];
         neural_all = [];
  
  
         % Flatten vertically along columns
         for i = 1:numel(drop_tst)
             dropsonde_all = [dropsonde_all; cat(2, drop_tst{i}(:))];
             lidar_all = [lidar_all; cat(2, lidar_tst{i}(:))];
             def_all = [def_all; cat(2, p2_tst{i}(:))]; 
             ensemble_all = [ensemble_all; cat(2, mlh_ens{i}(:))];
             linear_all = [linear_all; cat(2, mlh_linear{i}(:))];
             neural_all = [neural_all; cat(2, mlh_neural{i}(:))];
         end
  
         % Nans
         drop_nans = ~isnan(dropsonde_all);
         lidar_nans = ~isnan(lidar_all); 
         def_nans = ~isnan(def_all);
         ensemble_nans = ~isnan(ensemble_all);
         linear_nans = ~isnan(linear_all);
         neural_nans = ~isnan(neural_all);
         common_nans = drop_nans & lidar_nans & def_nans & ensemble_nans & linear_nans & neural_nans;

         %Lidar
         [nmb_lidar,nme_lidar,mb_lidar,me_lidar,rmse_lidar,corr_lidar] = calcStats_all_updated(common_nans,dropsonde_all,lidar_all); %default
  
         %Default 
         [nmb_def,nme_def,mb_def,me_def,rmse_def,corr_def] = calcStats_all_updated(common_nans,dropsonde_all,def_all);
  
         %Ensemble Tree Model
         [nmb_ens,nme_ens,mb_ens,me_ens,rmse_ens,corr_ens] = calcStats_all_updated(common_nans,dropsonde_all,ensemble_all);
  
         %Linear Model
         [nmb_linear,nme_linear,mb_linear,me_linear,rmse_linear,corr_linear] = calcStats_all_updated(common_nans,dropsonde_all,linear_all);
  
         %Neural Model 
         [nmb_neural,nme_neural,mb_neural,me_neural,rmse_neural,corr_neural] = calcStats_all_updated(common_nans,dropsonde_all,neural_all);
  
         %--Output to .txt file--
         NMB = [nmb_lidar, nmb_def, nmb_ens, nmb_linear, nmb_neural];
         NME = [nme_lidar, nme_def, nme_ens, nme_linear, nme_neural];
         MB = [mb_lidar, mb_def, mb_ens, mb_linear, mb_neural];
         ME = [me_lidar, me_def, me_ens, me_linear, me_neural];
         RMSE = [rmse_lidar, rmse_def, rmse_ens, rmse_linear, rmse_neural];
         CORR = [corr_lidar, corr_def, corr_ens, corr_linear, corr_neural];
  
  
         % Round the values to three decimal places
         NMB = round(NMB, 2);
         NME = round(NME, 2);
         MB = round(MB, 2);
         ME = round(ME, 2);
         RMSE = round(RMSE, 2);
         CORR = round(CORR, 2);
  
         columnNames = {'NMB', 'NME', 'MB', 'ME', 'RMSE', 'CORR'};
         rowNames = {'Lidar', 'Default', 'Ensemble', 'Linear', 'Neural'};
         tbl = table(NMB', NME', MB', ME', RMSE', CORR', 'VariableNames', columnNames, 'RowNames', rowNames);
  
  
  
         filename = 'DropStats_manuscript.txt';
         writetable(tbl, filename, 'Delimiter', '\t');
  
     end 
 
 end 
 












 
     
