text:
#heatmap plots for grouped and ungrouped genotype and read depth viz

"""
    genotype_heatmap2_new_legend(input::Array{Any,2},title::AbstractString,chrom_label_info,sample_names,chr_pos_tuple_list_rev,y_axis_label_option,x_axis_label_option)
generate heatmap of genotype data.
"""
function genotype_heatmap2_new_legend(input,title,chrom_label_info,sample_names,chr_pos_tuple_list_rev,y_axis_label_option,x_axis_label_option) #chr_pos_tuple_list_rev is rev because heatmap in plotly mirrors list for some reason.

            chr_pos_tuple_indices,chr_pos_tuple_list,sample_name_indices,sample_names,chrom_labels,chrom_label_indices,font_size = process_plot_inputs(chrom_label_info,sample_names,chr_pos_tuple_list_rev)
            title_no_underscores=replace(title, "_"=>' ')

             hover_text_array=generate_hover_text_array(chr_pos_tuple_list,sample_names,input,"GT")

             position_1,position_2,position_3,position_4,position_5,position_6,position_7,position_8,position_9,position_10,text_point_increment,legend_xaxis_anchor,legend_xaxis_anchor_2,x_axis_text_anchor,extra_right_margin_space,extra_right_margin_space_dp = generate_legend_increments_ungrouped(input)

               trace=heatmap(
                   z = input, x = 1:size(input, 2),y = 1:size(input, 1),

                   hoverinfo="text",
                   text=hover_text_array,

                   zauto=false,zmax=3,zmin=0,

                   transpose=true,
                   colorscale =    [[0, "rgb(255,255,255)"],
                                   [0.33, "rgb(51,106,145)"],
                                   [0.66, "rgb(65,165,137)"],
                                   [1, "rgb(251,231,65)"]],

                   gridcolor = "#E2E2E2",
                   showscale = false, #we define a custom legend using shapes and a text scatter trace
                   colorbar = attr(tickvals = [0, 1, 2, 3],
                   ticktext = ["No Call (0)", "Homozygous Reference", "Heterozygous Variant", "Homozygous Variant"])
                   );

                   #we define a custom legend using shapes and a text scatter trace #x_axis_text_anchor+extra_right_margin_space
                   trace2 = scatter(
                   ;x=[x_axis_text_anchor,x_axis_text_anchor,x_axis_text_anchor,x_axis_text_anchor,x_axis_text_anchor+extra_right_margin_space], y=y=[position_1+text_point_increment,position_3+text_point_increment,position_5+text_point_increment,position_7+text_point_increment,position_9+text_point_increment],
                     mode="text", name="legend label",
                     textposition="middle right",
                     text=["No Call", "Homozygous Reference", "Heterozygous Variant", "Homozygous Variant",""],
                     marker_size=300, textfont_family="Raleway, sans-serif")

                   shapes = [
                   rect(x0=legend_xaxis_anchor, y0=position_7, x1=legend_xaxis_anchor_2, y1=position_8, opacity=1.0, fillcolor="rgb(251,231,65)", line_color="black"),
                   rect(x0=legend_xaxis_anchor, y0=position_5, x1=legend_xaxis_anchor_2, y1=position_6, opacity=1.0, fillcolor="rgb(65,165,137)", line_color="black"),
                   rect(x0=legend_xaxis_anchor, y0= position_3, x1=legend_xaxis_anchor_2, y1=position_4, opacity=1.0, fillcolor="rgb(51,106,145)", line_color="black"),
                   rect(x0=legend_xaxis_anchor, y0= position_1, x1=legend_xaxis_anchor_2, y1=position_2, opacity=1.0, fillcolor="rgb(255,255,255)", line_color="black")]

                   data = [trace, trace2]

                if y_axis_label_option == "chromosomes"

                   layout = Layout(
                                   shapes=shapes,

                                   title = "$title_no_underscores",
                                   xaxis=attr(title="Sample ID", showgrid=false, zeroline=false, tickvals=sample_name_indices,
                                   ticktext=sample_names, tickangle=45,showticklabels=x_axis_label_option),
                                   yaxis=attr(title="Genomic Location", zeroline=false, tickvals=chrom_label_indices,
                                   ticktext=chrom_labels,tickfont_size=font_size,hovermode=true,automargin=true,showgrid=false)
                   )

                      plot(data,layout)

               elseif y_axis_label_option == "positions"

                   layout = Layout(shapes=shapes,
                                   title = "$title_no_underscores",
                                   xaxis=attr(title="Sample ID", showgrid=false, zeroline=false, tickvals=sample_name_indices,
                                   ticktext=sample_names, tickangle=45, showticklabels=x_axis_label_option),
                                   yaxis=attr(title="Genomic Location", zeroline=false, tickvals=chr_pos_tuple_indices,
                                   ticktext=chr_pos_tuple_list,tickfont_size=font_size,hovermode=true,automargin=true,showgrid=false)
                   )

                   plot(data,layout)

               elseif y_axis_label_option == "hover_positions"

                   layout = Layout(shapes=shapes,
                                   title = "$title_no_underscores",
                                   xaxis=attr(title="Sample ID", showgrid=false, zeroline=false, tickvals=sample_name_indices,
                                   ticktext=sample_names, tickangle=45,showticklabels=x_axis_label_option),
                                   yaxis=attr(title="Genomic Location", zeroline=false, tickvals=chr_pos_tuple_indices,
                                   ticktext=chr_pos_tuple_list,tickfont_size=font_size,hovermode=true,automargin=true,showticklabels=false)
                   )

                   plot(data,layout)

                else
                    println("--y_axis_labels is not a valid option. Choose positions or chromosomes")
                end
end

"""
   genotype_heatmap_with_groups(input::Array{Int64,2},title::String,chrom_label_info::Tuple{Array{String,1},Array{Int64,1},String},group_label_pack::Array{Any,1},id_list,chr_pos_tuple_list_rev,y_axis_label_option,trait_label_array,x_axis_label_option,number_rows)
generate heatmap of genotype data.
"""
function genotype_heatmap_with_groups(input::Array{Int64,2},title::String,chrom_label_info::Tuple{Array{String,1},Array{Int64,1},String},group_label_pack::Array{Any,1},id_list,chr_pos_tuple_list_rev,y_axis_label_option,trait_label_array,x_axis_label_option,number_rows)

    sample_name_indices,id_list,chrom_labels,chrom_label_indices,font_size,group_dividing_line,group1_label,group2_label,chr_pos_tuple_indices,chr_pos_tuple_list,font_size = process_plot_inputs_for_grouped_data(chrom_label_info::Tuple{Array{String,1},Array{Int64,1},String},group_label_pack::Array{Any,1},id_list,chr_pos_tuple_list_rev,trait_label_array)

    position_1,position_2,position_3,position_4,position_5,position_6,position_7,position_8,position_9,text_point_increment,trait_positions,legend_xaxis_anchor,legend_xaxis_anchor_2,x_axis_text_anchor,extra_right_margin_space = generate_legend_increments_grouped(input)

    position_1_trait_1,position_2_trait_1,position_1_trait_2,position_2_trait_2 = trait_positions

    title_no_underscores=replace(title, "_"=>' ')

    hover_text_array=generate_hover_text_array_grouped(chr_pos_tuple_list,id_list,input,"GT",number_rows)

                  trace = heatmap(
                       z = input, x=1:size(input, 2),y=1:size(input, 1),

                       hoverinfo="text",
                       text=hover_text_array,

                       zauto=false,zmax=3,zmin=-2,

                       transpose=true,

                                    colorscale = [
                                                  [0, "rgb(208, 211, 212)"], #light grey
                                                  [0.2, "rgb(151, 154, 154)"], #dark grey
                                                  [0.4, "rgb(255,255,255)"], #white
                                                  #[0.4, "rgb(56,25,90)"], #dark blue
                                                  [0.6, "rgb(51,106,145)"], #blue
                                                  [0.8, "rgb(65,165,137)"], #green
                                                  [1, "rgb(251,231,65)"] #yellow
                                                  ],

                      gridcolor = "#E2E2E2",
                      showscale = false, #we define a custom legend using shapes and a text scatter trace
                      colorbar = attr(tickvals = [-2, -1, 0, 1, 2, 3],
                      ticktext = ["Trait 2", "Trait 1", "No Call", "Homozygous Reference", "Heterozygous Variant", "Homozygous Variant"])
                      );

                    #we define a custom legend using shapes and a text scatter trace

                    trace2 = scatter(
                    ;x=[x_axis_text_anchor,x_axis_text_anchor,x_axis_text_anchor,x_axis_text_anchor,x_axis_text_anchor+extra_right_margin_space,x_axis_text_anchor,x_axis_text_anchor], y=y=[position_1+text_point_increment,position_2+text_point_increment,position_3+text_point_increment,position_4+text_point_increment,position_5+text_point_increment,position_7+text_point_increment,position_8+text_point_increment],
                      mode="text", name="legend label",
                      textposition="middle right",
                      text=["No Call", "Homozygous Reference", "Heterozygous Variant", "Homozygous Variant"," ", "Trait 2","Trait 1"],
                      marker_size=300, textfont_family="Raleway, sans-serif")

                    shapes = [
                    rect(x0=legend_xaxis_anchor, y0=position_4, x1=legend_xaxis_anchor_2, y1=position_5, opacity=1.0, fillcolor="rgb(251,231,65)", line_color="black"),
                    rect(x0=legend_xaxis_anchor, y0=position_3, x1=legend_xaxis_anchor_2, y1=position_4, opacity=1.0, fillcolor="rgb(65,165,137)", line_color="black"),
                    rect(x0=legend_xaxis_anchor, y0= position_2, x1=legend_xaxis_anchor_2, y1=position_3, opacity=1.0, fillcolor="rgb(51,106,145)", line_color="black"),
                    rect(x0=legend_xaxis_anchor, y0= position_1, x1=legend_xaxis_anchor_2, y1=position_2, opacity=1.0, fillcolor="rgb(255,255,255)", line_color="black"),
                    rect(x0=legend_xaxis_anchor, y0= position_7, x1=legend_xaxis_anchor_2, y1=position_8, opacity=1.0, fillcolor="rgb(208, 211, 212)", line_color="black"),
                    rect(x0=legend_xaxis_anchor, y0= position_8, x1=legend_xaxis_anchor_2, y1=position_9, opacity=1.0, fillcolor="rgb(151, 154, 154)", line_color="black")
                    ]

                     data = [trace, trace2]

      if y_axis_label_option == "chromosomes"

          layout = Layout(
                          title = "$title_no_underscores",

                          xaxis=attr(
                          title="Sample ID (Grouped by $(group1_label) | $(group2_label))",
                          showgrid=false,
                          zeroline=false,
                          tickvals=sample_name_indices,
                          ticktext=id_list,
                          tickangle=45,
                          showticklabels=x_axis_label_option),

                          yaxis=attr(
                          title="Genomic Location",
                          zeroline=false,
                          tickvals=chrom_label_indices,
                          ticktext=chrom_labels,
                          tickfont_size=font_size,
                          hovermode=true,
                          automargin=true,showgrid=false),

                          shapes=shapes,
                          yaxis_range = [1:size(input,1)],
                          xaxis_range = [1:size(input,2)]
                        )

                        plot(data,layout)

      elseif y_axis_label_option == "positions"

          layout = Layout(
                          title = "$title_no_underscores",

                          xaxis=attr(title="Sample ID (Grouped by $(group1_label) | $(group2_label))", showgrid=false, zeroline=false, tickvals=sample_name_indices,
                          ticktext=id_list, tickangle=45,showticklabels=x_axis_label_option),

                          yaxis=attr(title="Genomic Location", zeroline=false, tickvals=chr_pos_tuple_indices,
                          ticktext=chr_pos_tuple_list,tickfont_size=font_size,hovermode=true,automargin=true,showgrid=false),

                          shapes=shapes, yaxis_range = [1:size(input,1)], xaxis_range = [1:size(input,2)] #line bright white
                        )

                        plot(data,layout)

      elseif y_axis_label_option == "hover_positions"

        layout = Layout(
                        title = "$title_no_underscores",
                        xaxis=attr(title="Sample ID (Grouped by $(group1_label) | $(group2_label))", showgrid=false, zeroline=false, tickvals=sample_name_indices,
                        ticktext=id_list, tickangle=45,showticklabels=x_axis_label_option),

                        yaxis=attr(title="Genomic Location", zeroline=false, tickvals=chr_pos_tuple_indices,
                        ticktext=chr_pos_tuple_list,tickfont_size=font_size,hovermode=true,automargin=true,showticklabels=false),

                        shapes=shapes, yaxis_range = [1:size(input,1)], xaxis_range = [1:size(input,2)]
        )

                        plot(data,layout)

      else
          println("--y_axis_labels is not a valid option. Choose positions or chromosomes")
      end
  end

"""
    dp_heatmap2(input::Array{Int64,2},title::String,chrom_label_info::Tuple{Array{String,1},Array{Int64,1},String}, sample_names,chr_pos_tuple_list_rev,y_axis_label_option,x_axis_label_option)
generate heatmap of read depth data.
"""
function dp_heatmap2(input, title, chrom_label_info, sample_names,chr_pos_tuple_list_rev,y_axis_label_option,x_axis_label_option)

    chr_pos_tuple_indices,chr_pos_tuple_list,sample_name_indices,sample_names,chrom_labels,chrom_label_indices,font_size = process_plot_inputs(chrom_label_info,sample_names,chr_pos_tuple_list_rev)

    title_no_underscores=replace(title, "_"=>' ')

    hover_text_array=generate_hover_text_array(chr_pos_tuple_list,sample_names,input,"DP")

    position_1,position_2,position_3,position_4,position_5,position_6,position_7,position_8,position_9,position_10,text_point_increment,legend_xaxis_anchor,legend_xaxis_anchor_2,x_axis_text_anchor,extra_right_margin_space,extra_right_margin_space_dp = generate_legend_increments_ungrouped(input)

    trace=heatmap(
        z = input, x=1:size(input, 2),y=1:size(input, 1),

       hoverinfo="text",
       text=hover_text_array,

        zauto=false,zmax=100,zmin=-20,

        transpose=true,
        colorscale = [
                     [0, "rgb(255,255,255)"],
                     [0.2, "rgb(153,231,255)"],
                     [0.3, "rgb(79,146,255)"],
                     [0.41, "rgb(43,124,255)"],
                     [1, "rgb(0,64,168)"]
                     ],

        colorbar = attr(tickvals = [-20,0,20,40,60,80,99],
        title="Read Depth",
        ticktext = ["No Call","0","20","40","60","80","100+"]),

        gridcolor = "#E2E2E2",
        showscale = false, #we define a custom legend using shapes and a text scatter trace
        );

        #we define a custom legend using shapes and a text scatter trace
        trace2 = scatter(
        ;x=[x_axis_text_anchor,x_axis_text_anchor,x_axis_text_anchor,x_axis_text_anchor,x_axis_text_anchor,x_axis_text_anchor+extra_right_margin_space_dp],
        y=[position_1+text_point_increment,position_2+text_point_increment,position_3+text_point_increment,position_4+text_point_increment,position_5+text_point_increment,position_6+text_point_increment],
          mode="text", name="legend label",
          textposition="middle right",
          text=["0", "10", "30", "70","100+",""],
          marker_size=300, textfont_family="Raleway, sans-serif")

        shapes = [
        rect(x0=legend_xaxis_anchor, y0=position_5, x1=legend_xaxis_anchor_2, y1=position_6, opacity=1.0, fillcolor="rgb(0,64,168)", line_color="black"),
        rect(x0=legend_xaxis_anchor, y0=position_4, x1=legend_xaxis_anchor_2, y1=position_5, opacity=1.0, fillcolor="rgb(43,124,255)", line_color="black"),
        rect(x0=legend_xaxis_anchor, y0=position_3, x1=legend_xaxis_anchor_2, y1=position_4, opacity=1.0, fillcolor="rgb(79,146,255)", line_color="black"),
        rect(x0=legend_xaxis_anchor, y0= position_2, x1=legend_xaxis_anchor_2, y1=position_3, opacity=1.0, fillcolor="rgb(153,231,255)", line_color="black"),
        rect(x0=legend_xaxis_anchor, y0= position_1, x1=legend_xaxis_anchor_2, y1=position_2, opacity=1.0, fillcolor="rgb(255,255,255)", line_color="black")
        ]

        if y_axis_label_option == "chromosomes"

           layout = Layout(
                           shapes=shapes,
                           title = "$title_no_underscores",
                           xaxis=attr(title="Sample ID", showgrid=false, zeroline=false, tickvals=sample_name_indices,
                           ticktext=sample_names, tickangle=45,showticklabels=x_axis_label_option),
                           yaxis=attr(title="Genomic Location", zeroline=false, tickvals=chrom_label_indices,
                           ticktext=chrom_labels,tickfont_size=font_size,hovermode=true,automargin=true,showgrid=false)
           )

          data = [trace, trace2]
              plot(data,layout)

       elseif y_axis_label_option == "positions"

           layout = Layout(
                           shapes=shapes,
                           title = "$title_no_underscores",
                           xaxis=attr(title="Sample ID", showgrid=false, zeroline=false, tickvals=sample_name_indices,
                           ticktext=sample_names, tickangle=45,showticklabels=x_axis_label_option),
                           yaxis=attr(title="Genomic Location", zeroline=false, tickvals=chr_pos_tuple_indices,
                           ticktext=chr_pos_tuple_list,tickfont_size=font_size,hovermode=true,automargin=true,showgrid=false)
           )

      data = [trace, trace2]
           plot(data,layout)

       elseif y_axis_label_option == "hover_positions"

           layout = Layout(
                           shapes=shapes,
                           title = "$title_no_underscores",
                           xaxis=attr(title="Sample ID", showgrid=false, zeroline=false, tickvals=sample_name_indices,
                           ticktext=sample_names, tickangle=45,showticklabels=x_axis_label_option),
                           yaxis=attr(title="Genomic Location", zeroline=false, tickvals=chr_pos_tuple_indices,
                           ticktext=chr_pos_tuple_list,tickfont_size=font_size,hovermode=true,automargin=true,showticklabels=false)
           )
          data = [trace, trace2]
               plot(data,layout)

        else
            println("--y_axis_labels is not a valid option. Choose positions or chromosomes")
        end
end

"""
    dp_heatmap2_with_groups(input::Array{Int64,2},title::String,chrom_label_info::Tuple{Array{String,1},Array{Int64,1},String},group_label_pack::Array{Any,1},id_list,chr_pos_tuple_list_rev,y_axis_label_option,trait_label_array,x_axis_label_option,number_rows)
generate heatmap of read depth data with grouped samples.
"""
function dp_heatmap2_with_groups(input::Array{Int64,2},title::String,chrom_label_info::Tuple{Array{String,1},Array{Int64,1},String},group_label_pack::Array{Any,1},id_list,chr_pos_tuple_list_rev,y_axis_label_option,trait_label_array,x_axis_label_option,number_rows)

    sample_name_indices,id_list,chrom_labels,chrom_label_indices,font_size,group_dividing_line,group1_label,group2_label,chr_pos_tuple_indices,chr_pos_tuple_list,font_size = process_plot_inputs_for_grouped_data(chrom_label_info::Tuple{Array{String,1},Array{Int64,1},String},group_label_pack::Array{Any,1},id_list,chr_pos_tuple_list_rev,trait_label_array)

    title_no_underscores=replace(title, "_"=>' ')

    hover_text_array=generate_hover_text_array_grouped(chr_pos_tuple_list,id_list,input,"DP",number_rows)

    position_1,position_2,position_3,position_4,position_5,position_6,position_7,position_8,position_9,text_point_increment,trait_positions,legend_xaxis_anchor,legend_xaxis_anchor_2,x_axis_text_anchor,extra_right_margin_space = generate_legend_increments_grouped(input)

    position_1_trait_1,position_2_trait_1,position_1_trait_2,position_2_trait_2 = trait_positions

    trace=heatmap(

        z = input, x=1:size(input, 2),y=1:size(input, 1),

        hoverinfo="text",
        text=hover_text_array,

        zauto=false,zmax=100,zmin=-60,

        transpose=true,

                     colorscale = [
                                      [0, "rgb(208, 211, 212)"], #light grey
                                      [0.125, "rgb(151, 154, 154)"], #dark grey
                                      [0.25, "rgb(255, 255, 255)"],
                                      [0.4, "rgb(153,231,255)"],
                                      [0.5, "rgb(79,146,255)"],
                                      [0.5625, "rgb(43,124,255)"],
                                      [1, "rgb(0,64,168)"]
                                  ],

                     colorbar = attr(tickvals = [-60, -40, -20, 0, 20, 40, 60, 80, 99],
                     title="Depth / Trait",
                     ticktext = ["Trait 2","Trait 1","No Call","0","20","40","60","80","100+"]),
                     gridcolor = "#E2E2E2",
                     showscale = false, #we define a custom legend using shapes and a text scatter trace
                     );

         #we define a custom legend using shapes and a text scatter trace
         trace2 = scatter(
         ;x=[x_axis_text_anchor,x_axis_text_anchor,x_axis_text_anchor,x_axis_text_anchor,x_axis_text_anchor,x_axis_text_anchor+extra_right_margin_space,x_axis_text_anchor,x_axis_text_anchor],
         y=[position_1+text_point_increment,position_2+text_point_increment,position_3+text_point_increment,position_4+text_point_increment,position_5+text_point_increment,position_6 + text_point_increment, position_7+text_point_increment,position_8+text_point_increment],
           mode="text", name="legend label",
           textposition="middle right",
           text=["No Call", "0", "10", "20","100+","","Trait 2","Trait 1"],
           marker_size=300, textfont_family="Raleway, sans-serif")

         data = [trace, trace2]

         shapes = [
         rect(x0=legend_xaxis_anchor, y0= position_7, x1=legend_xaxis_anchor_2, y1=position_8, opacity=1.0, fillcolor="rgb(208, 211, 212)", line_color="black"),
         rect(x0=legend_xaxis_anchor, y0= position_8, x1=legend_xaxis_anchor_2, y1=position_9, opacity=1.0, fillcolor="rgb(151, 154, 154)", line_color="black"),
         rect(x0=legend_xaxis_anchor, y0=position_5, x1=legend_xaxis_anchor_2, y1=position_6, opacity=1.0, fillcolor="rgb(37,114,242)", line_color="black"),
         rect(x0=legend_xaxis_anchor, y0=position_4, x1=legend_xaxis_anchor_2, y1=position_5, opacity=1.0, fillcolor="rgb(67,138,254)", line_color="black"),
         rect(x0=legend_xaxis_anchor, y0=position_3, x1=legend_xaxis_anchor_2, y1=position_4, opacity=1.0, fillcolor="rgb(110,182,255)", line_color="black"),
         rect(x0=legend_xaxis_anchor, y0= position_2, x1=legend_xaxis_anchor_2, y1=position_3, opacity=1.0, fillcolor="rgb(153,231,255)", line_color="black"),
         rect(x0=legend_xaxis_anchor, y0= position_1, x1=legend_xaxis_anchor_2, y1=position_2, opacity=1.0, fillcolor="rgb(255,255,255)", line_color="black")
         ]

if y_axis_label_option == "chromosomes"

    layout = Layout(
                    shapes=shapes,
                    title = "$title_no_underscores",

                    xaxis=attr(title="Sample ID (Grouped by $(group1_label) | $(group2_label))",
                    showgrid=false, zeroline=false, tickvals=sample_name_indices,
                    ticktext=id_list, tickangle=45,showticklabels=x_axis_label_option),

                    yaxis=attr(title="Genomic Location", zeroline=false,
                    tickvals=chrom_label_indices,ticktext=chrom_labels,
                    tickfont_size=font_size,hovermode=true,automargin=true,showgrid=false),

                    yaxis_range = [1:size(input,1)],
                    xaxis_range = [1:size(input,2)]
                  )

      plot(data,layout)

elseif y_axis_label_option == "positions"

    layout = Layout(
                    shapes=shapes,
                    title = "$title_no_underscores",

                    xaxis=attr(title="Sample ID (Grouped by $(group1_label) | $(group2_label))",
                    showgrid=false, zeroline=false, tickvals=sample_name_indices,
                    ticktext=id_list, tickangle=45,showticklabels=x_axis_label_option),

                    yaxis=attr(title="Genomic Location", zeroline=false,
                    tickvals=chr_pos_tuple_indices,
                    ticktext=chr_pos_tuple_list,tickfont_size=font_size,
                    hovermode=true,automargin=true,showgrid=false),

                    yaxis_range = [1:size(input,1)],
                    xaxis_range = [1:size(input,2)]
                  )

      plot(data,layout)

elseif y_axis_label_option == "hover_positions"

  layout = Layout(
                  shapes=shapes,
                  title = "$title_no_underscores",
                  xaxis=attr(title="Sample ID (Grouped by $(group1_label) | $(group2_label))", showgrid=false, zeroline=false, tickvals=sample_name_indices,
                  ticktext=id_list, tickangle=45,showticklabels=x_axis_label_option),

                  yaxis=attr(title="Genomic Location", zeroline=false, tickvals=chr_pos_tuple_indices,
                  ticktext=chr_pos_tuple_list,tickfont_size=font_size,hovermode=true,automargin=true,showticklabels=false),

                  yaxis_range = [1:size(input,1)], xaxis_range = [1:size(input,2)] #line bright white
  )
  plot(data,layout)

else
    println("--y_axis_labels is not a valid option. Choose positions or chromosomes")
end

end

#Scatter plots for average read depth viz

"""
    avg_sample_dp_scatter(sample_avg_list::Array{Float64,1},sample_names,x_axis_label_option)
generate line chart of average read depths of each sample.
"""
function avg_sample_dp_scatter(sample_avg_list::Array{Float64,1},sample_names,x_axis_label_option)

    sample_name_indices = collect(1:1:size(sample_names,2))

    trace = scatter(;x=1:size(sample_avg_list,1), y=sample_avg_list,mode="markers")

    layout = Layout(title="Average Sample Read Depth",
    xaxis=attr(title="Samples",ticktext=sample_names,tickvals=sample_name_indices,
    tickangle=45,tickfont_size=5,showticklabels=x_axis_label_option),
    yaxis=attr(title="Average Read Depth"),showticklabels=false)

    plot(trace,layout)
end

"""
    avg_variant_dp_line_chart(variant_avg_list::Array{Float64,1},chr_pos_tuple_list,y_axis_label_option,chrom_label_info)
generate line chart of average read depths of each variant.
"""
function avg_variant_dp_line_chart(variant_avg_list::Array{Float64,1},chr_pos_tuple_list,y_axis_label_option,chrom_label_info)

    chr_pos_tuple_indices = collect(1:1:size(chr_pos_tuple_list,1))
    chrom_labels = chrom_label_info[1]
    returnXY_column1!(chrom_labels)
    chrom_label_indices = chrom_label_info[2]
    font_size = chrom_label_info[3]

    trace = scatter(;x=1:size(variant_avg_list,1), y=variant_avg_list,mode="markers")

    if y_axis_label_option == "positions"
        layout = Layout(title="Average Variant Read Depth",xaxis=attr(title="Variant Positions",
        ticktext=chr_pos_tuple_list,
        tickvals=chr_pos_tuple_indices,
        tickangle=45,showticklabels=true),yaxis=attr(title="Average Read Depth"))
    elseif y_axis_label_option == "chromosomes"
        layout = Layout(title="Average Variant Read Depth",xaxis=attr(title="Variant Positions",
        tickvals=chrom_label_indices,
        ticktext=chrom_labels,
        showticklabels=true),yaxis=attr(title="Average Read Depth"))
    elseif y_axis_label_option == "hover_positions"
        layout = Layout(title="Average Variant Read Depth",xaxis=attr(title="Variant Positions",
        ticktext=chr_pos_tuple_list,
        tickvals=chr_pos_tuple_indices,
        tickangle=45,showticklabels=false,hovermode=true),yaxis=attr(title="Average Read Depth"))
    else
        println("--y_axis_labels is not a valid option. Choose positions or chromosomes")
    end

    plot(trace,layout)

end

"""
    process_plot_inputs(chrom_label_info,sample_names,chr_pos_tuple_list_rev)
Prepares input for heatmap plot function for both genotype and read depth plots without --group_samples flag.
"""
function process_plot_inputs(chrom_label_info,sample_names,chr_pos_tuple_list_rev)

    chr_pos_tuple_indices = collect(1:1:size(chr_pos_tuple_list_rev,1))

    #chr_pos_tuple_list=Array{Tuple{Any,Int64}}(undef,0) #use this if tuple to string doesnt work #version
    chr_pos_tuple_list=Array{String}(undef,0)

    for i in chr_pos_tuple_list_rev
        push!(chr_pos_tuple_list, i)
    end

    sample_name_indices = collect(1:1:size(sample_names,2))

    chrom_labels = chrom_label_info[1]
    returnXY_column1!(chrom_labels)
    chrom_label_indices = chrom_label_info[2]
    font_size = chrom_label_info[3]

    return chr_pos_tuple_indices,chr_pos_tuple_list,sample_name_indices,sample_names,chrom_labels,chrom_label_indices,font_size
end

"""
    process_plot_inputs_for_grouped_data(chrom_label_info::Tuple{Array{String,1},Array{Int64,1},String},group_label_pack::Array{Any,1},id_list,chr_pos_tuple_list_rev,trait_label_array)
Prepares input for heatmap plot function for both genotype and read depth plots with --group_samples flag.
"""
function process_plot_inputs_for_grouped_data(chrom_label_info::Tuple{Array{String,1},Array{Int64,1},String},group_label_pack::Array{Any,1},id_list,chr_pos_tuple_list_rev,trait_label_array)

    sample_name_indices = collect(1:1:size(id_list,2))

    chrom_labels = chrom_label_info[1]
    returnXY_column1!(chrom_labels)
    chrom_label_indices = chrom_label_info[2]
    font_size = chrom_label_info[3]

    group_dividing_line=group_label_pack[3]
    group1_label=group_label_pack[4]
    group2_label=group_label_pack[5]

    chr_pos_tuple_list=Array{String}(undef,0)

    for i in chr_pos_tuple_list_rev
        push!(chr_pos_tuple_list, i)
    end

    y_labels=vcat(chr_pos_tuple_list,trait_label_array)

    y_labels_indices = collect(1:1:size(y_labels,1))

    return sample_name_indices,id_list,chrom_labels,chrom_label_indices,font_size,group_dividing_line,group1_label,group2_label,y_labels_indices,y_labels,font_size

end

"""
     generate_legend_increments_ungrouped(input)
Dynamically generates positons for shapes that build categorical colorscale.
"""
function generate_legend_increments_ungrouped(input)

    number_rows = size(input,1)
    increment = number_rows/30

    number_samples = size(input,2)
    distance_from_plot_edge = number_samples/40
    width_of_color_box = number_samples/40
    legend_xaxis_anchor = number_samples + distance_from_plot_edge
    legend_xaxis_anchor_2 = legend_xaxis_anchor + width_of_color_box
    x_axis_text_anchor = legend_xaxis_anchor_2 + number_samples/40
    extra_right_margin_space=number_samples/4
    extra_right_margin_space_dp = number_samples/10

    position_1 = number_rows/2
    position_2 = position_1 + increment
    position_3 = position_2 + increment
    position_4 = position_3 + increment
    position_5 = position_4 + increment
    position_6 = position_5 + increment
    position_7 = position_6 + increment
    position_8 = position_7 + increment
    position_9 = position_8 + increment
    position_10 = position_9 + increment
    text_point_increment = increment/2
    #trait_positions = position_1_trait_1,position_2_trait_1,position_1_trait_2,position_1_trait_2

    return position_1,position_2,position_3,position_4,position_5,position_6,position_7,position_8,position_9,position_10,text_point_increment,legend_xaxis_anchor,legend_xaxis_anchor_2,x_axis_text_anchor,extra_right_margin_space,extra_right_margin_space_dp
end

"""
     generate_legend_increments_grouped(input)
Dynamically generates positons for shapes that build categorical colorscale including two color boxes for traits 1 and 2
"""
function generate_legend_increments_grouped(input)

    number_rows = size(input,1)
    increment = number_rows/30

    number_samples = size(input,2)
    distance_from_plot_edge = number_samples/40
    width_of_color_box = number_samples/40
    legend_xaxis_anchor = number_samples + distance_from_plot_edge
    legend_xaxis_anchor_2 = legend_xaxis_anchor + width_of_color_box
    x_axis_text_anchor = legend_xaxis_anchor_2 + number_samples/40
    extra_right_margin_space = number_samples/3

    position_1 = number_rows/3
    position_2 = position_1 + increment
    position_3 = position_2 + increment
    position_4 = position_3 + increment
    position_5 = position_4 + increment
    position_6 = position_5 + increment
    position_7 = position_6 + increment
    position_8 = position_7 + increment
    position_9 = position_8 + increment
    #position_10 = position_9 + increment
    text_point_increment = increment/2
    position_1_trait_1 = position_9 + increment
    position_2_trait_1 = position_1_trait_1 + increment
    position_1_trait_2 = position_2_trait_1 + increment
    position_2_trait_2 = position_1_trait_2 + increment
    trait_positions = position_1_trait_1,position_2_trait_1,position_1_trait_2,position_2_trait_2

    return position_1,position_2,position_3,position_4,position_5,position_6,position_7,position_8,position_9,text_point_increment,trait_positions,legend_xaxis_anchor,legend_xaxis_anchor_2,x_axis_text_anchor,extra_right_margin_space
end


text:
#=
This is a script that searches some file for a string and replaces
said string with another. In this case, the file is usually a
chemical network and the string is a species name.
=#

species = ARGS[1]
println("Species is $species")

newspecies = ARGS[2]
println("New species is $newspecies")

if ( ARGS[3] == "grain" )
  basefile = "grain_reactions.in"
  network = open("$(species)_$(basefile)")
  grainsp = "J$species"
  newgrain = "J$(newspecies)"
elseif ( ARGS[3] == "gas" )
  basefile = "gas_reactions.in"
  network = open("$(species)_$(basefile)")
end

outfile = open("$(newspecies)_$(basefile)","w")

for ln in eachline(network)
  newln = replace(ln,species,newspecies)
# apparently, replace() is smart enough to
# not care if there is a J in front of the
# string
#  if ( ARGS[3] == "grain" )
#    newln = replace(newln,grainsp,newgrain)
#  end
  println(newln)
  write(outfile,newln)
end

close(network)
close(outfile)

ext: jl


text:
# These examples illustrate common operations on polyhedra using [Polyhedra.jl](https://github.com/JuliaPolyhedra/Polyhedra.jl):
#
# - [Convex hull](#Convex-hull).
# - [Intersection](#Intersection).
#
# We start by choosing a polyhedral library that will be used for computing the H-representation from the V-representation and vice-versa as well as removing redundant points.
# In these example, we use the default library available in Polyhedra but it can be replaced by any other library listed [here](https://juliapolyhedra.github.io/), e.g. by changing the last two lines below by `import CDDLib` and `lib = CDDLib.Library()` to use [CDDLib](https://github.com/JuliaPolyhedra/CDDLib.jl).

using Test #jl
using Polyhedra
import GLPK
lib = DefaultLibrary{Float64}(GLPK.Optimizer)

# ### Convex hull
#
# The binary convex hull operation between two polyhedra is obtained with the [`convexhull`](https://juliapolyhedra.github.io/Polyhedra.jl/latest/utilities.html#Polyhedra.convexhull) function.
#
# Below we compute the convex hull of the union of two polygons from their V-representation.

P1 = polyhedron(vrep([
    -1.9 -1.7
    -1.8  0.5
     1.7  0.7
     1.9 -0.3
     0.9 -1.1
]), lib)

P2 = polyhedron(vrep([
    -2.5 -1.1
    -0.8  0.8
     0.1  0.9
     1.8 -1.2
     1.3  0.1
]), lib)

Pch = convexhull(P1, P2)

# Note that the convex hull operation is done in the V-representation so no representation conversion is needed for this operation since `P1` and `P2` where constructed from their V-representation:

hrepiscomputed(P1), hrepiscomputed(P2), hrepiscomputed(Pch) #!jl
@test (hrepiscomputed(P1), hrepiscomputed(P2), hrepiscomputed(Pch)) == (false, false, false) #jl

# Let us note that the `convexhull` of a V-representation contains points and rays and represents the convex hull of the points together with the conic hull of the rays. So, `convexhull(P1, P2)` does the union of the vertices:

npoints(Pch) #!jl
@test npoints(Pch) == 10 #jl

# However, if we want to remove the redundant points we can use `removevredundancy!`:

removevredundancy!(Pch)
npoints(Pch) #!jl
@test npoints(Pch) == 8 #jl

# We can plot the polygons and the convex hull of their union using the `plot` function. For further plotting options see the [Plots.jl](http://docs.juliaplots.org/latest/) documentation.
# We can see below the 8 redundant points highlighted with green dots, the two points that are not highlighted are the redundant ones.

using Plots                          #!jl
plot(P1, color="blue", alpha=0.2)    #!jl
plot!(P2, color="red", alpha=0.2)    #!jl
plot!(Pch, color="green", alpha=0.1) #!jl
scatter!(Pch, color="green")         #!jl

# ### Intersection
#
# Intersection of polyhedra is obtained with the [`intersect`](https://juliapolyhedra.github.io/Polyhedra.jl/latest/utilities.html#Base.intersect) function.
#
# Below we compute the intersection of the two polygons from the previous example.

Pint = intersect(P1, P2)
@test nhalfspaces(Pint) == 10 #jl

# While `P1` and `P2` have been constructed from their V-representation, their H-representation has been computed to build the intersection `Pint`.

hrepiscomputed(P1), vrepiscomputed(P1), hrepiscomputed(P2), vrepiscomputed(P2) #!jl
@test (hrepiscomputed(P1), vrepiscomputed(P1), hrepiscomputed(P2), vrepiscomputed(P2)) == (true, true, true, true) #jl

# On the other hand, `Pint` is constructed from its H-representation hence its V-representation has not been computed yet.

hrepiscomputed(Pint), vrepiscomputed(Pint) #!jl
@test (hrepiscomputed(Pint), vrepiscomputed(Pint)) == (true, false) #jl

# We can obtain the number of points in the intersection with `npoints` as follows:

npoints(Pint) #!jl
@test npoints(Pint) == 8 #jl

# Note that this triggers the computation of the V-representation:

hrepiscomputed(Pint), vrepiscomputed(Pint) #!jl
@test (hrepiscomputed(Pint), vrepiscomputed(Pint)) == (true, true) #jl

# We can plot the polygons and their intersection using the `plot` function.

using Plots                            #!jl
plot(P1, color="blue", alpha=0.2)      #!jl
plot!(P2, color="red", alpha=0.2)      #!jl
plot!(Pint, color="yellow", alpha=0.6) #!jl

max_stars_repo_path: examples/Convex hull and intersection.jl


text:
module GPUArrays

export GPUArray, gpu_call, thread_blocks_heuristic, global_size, synchronize_threads
export linear_index, @linearidx, @cartesianidx, convolution!, device, synchronize
export JLArray

using Serialization
using Random
using LinearAlgebra
using Printf

using LinearAlgebra.BLAS
using Base.Cartesian

using FFTW

using Adapt

include("abstractarray.jl")
include("abstract_gpu_interface.jl")
include("ondevice.jl")
include("base.jl")
include("construction.jl")
include("blas.jl")
include("broadcast.jl")
include("devices.jl")
include("heuristics.jl")
include("indexing.jl")
include("linalg.jl")
include("mapreduce.jl")
include("vectors.jl")
include("convolution.jl")
include("random.jl")

include("array.jl")

include("testsuite.jl")

end # module


text:
#= ercot_stoch_reliability.jl
solves the a stochastic two stage unit-commitment model in which slow generators
are committed in the first stage and fast generators can be committed in real time.
Uses formulation based on Papavasiliou and Oren (2013).
Demand response can function as a slow generator with an advance commmitment
and an advane production schedule if desired.
Uncertainty comes from total demand.

cmd line format should be:
julia ercot_stoch_reliability.jl <date> <inputs_file_name.csv> <multi-runTF> <period_name>

multi-runTF is usually true: it indicates that the current simulation is one time
period of a group of periods that make up a single simulation. This just affects
the file structure.

Written under Julia 0.6.4
Patricia Levi
pjlevi@stanford.edu
NOV 2018
=#

### Packages ###
using JuMP
using Gurobi
using DataFrames
using CSV
# using JLD2
include("convert3dto2d.jl")
include("make_scenarios.jl")
include("writecsvmulti.jl")

test = Gurobi.Env() # test that gurobi is working

# -------------------------------------------
# USER CONTROLS & CMD LINE ARGS ##########
# -------------------------------------------

# USER PARAMS #
no_vars = false #stops execution before making variables
debug = false  # stops execution before solving model

sherlock_fol = "~/dr_stoch_uc/"
sherlock_input_file = "inputs/"
sherlock_output_file = "outputs/"
sherlock_params_fol = "code/"
laptop_fol = "~/Documents/dr_stoch_uc/"
laptop_input_file = "inputs/"
laptop_output_file = "outputs/"
laptop_params_fol = "code/"

dr_gens = ["DR"]
notdr_gens = ["BIOMASS","COAL","GAS","GAS_CC","GAS_CT","GAS_ICE","GAS_ST","HYDRO",
            "LANDFILL_GAS","NUCLEAR","OIL","SOLAR","WIND"    ]

####### END USER CONTROLS ##########

# --------------------------------------------------------------------------------------
# FILEPATH SETUP
# -------------------------------------------
if split(pwd(),"/")[2] == "Users"
    Sherlock = false
else
    Sherlock = true # on sherlock? where are folders?
end
@show Sherlock

if Sherlock
    base_fol = sherlock_fol
    input_fol = string(sherlock_fol,sherlock_input_file)
    params_fol = string(sherlock_fol, sherlock_params_fol)
else
    base_fol = laptop_fol
    input_fol = string(laptop_fol,laptop_input_file)
    params_fol = string(laptop_fol, laptop_params_fol)
end
default_data_fol = string(input_fol,"ercot_default/")

# PARSE CMD LINE ARGS #
ARGNAMES = ["date" ,"inputs_file_name","input_verion" ,"multi-runTF", "period_name" ]
defaultARGS = [Dates.format(Dates.now(),"Y-m-d"),"inputs_ercot.csv","tiny","true","periods_1_5353_5400.csv"]
localARGS = length(ARGS) > 0 ? ARGS : defaultARGS #if ARGS supplied, use those. otherwise, use default
nargs = length(localARGS)
@show localARGS

if nargs == 5
    submitdate = localARGS[1]
    input_file_name = localARGS[2]
    input_version = localARGS[3]
    multiTF = parse(Bool,lowercase(localARGS[4]))
    periodID = localARGS[5]
elseif nargs > 5
    error(string("Too many arguments supplied. Need ",join(ARGNAMES[1,:]," ")))
elseif nargs <5
    warn("not enough arguments supplied. Need ", join(ARGNAMES[1,:]," "))
end

read_inputs = CSV.read(string(params_fol,input_file_name))
read_inputs = read_inputs[:,[:input_name, Symbol(input_version)]]
# add cmd line args to read_inputs file
newdf = DataFrame(input_name = ARGNAMES[1:length(localARGS)], args = localARGS)
names!(newdf.colindex,map(parse,["input_name",input_version]))
read_inputs = vcat(read_inputs, newdf)
@show read_inputs
# ~ transpose read_inputs so that values can be referenced by name
read_inputs[:,:rowkey] = 1
inputs = unstack(read_inputs,:rowkey, :input_name, Symbol(input_version))

# parse out non-string inputs
startlim = parse(Float64,inputs[1,:startlim])
hourlim = parse(Float64,inputs[1,:hourlim])
energylim = parse(Float64,inputs[1,:energylim])
ramplims = parse(Float64,inputs[1,:ramplims])
durationlim = parse(Float64, inputs[1,:durationlim])
dr_override = parse(Bool,lowercase(inputs[1,:dr_override]))
dr_varcost = parse(Float64,inputs[1,:dr_varcost])
randScenarioSel = parse(Bool,lowercase(inputs[1,:randScenarioSel]))
trueBinaryStartup = parse(Bool,lowercase(inputs[1,:trueBinaryStartup]))
DRtype = parse(Int64,inputs[1,:DRtype])
nrandp = parse(Int64,inputs[1,:nrandp])
int_length = parse(Int64,inputs[1,:intlength])
debug = !parse(Bool, lowercase(inputs[1,:solve_model]))
MIPFocusParam = parse(Int64,inputs[1,:MIPFocusParam])
MIPGapParam = parse(Float64,inputs[1,:MIPGapParam])
NodefileStartParam = parse(Float64,inputs[1,:NodefileStartParam])
MethodParam = parse(Int64,inputs[1,:MethodParam])

if Sherlock
    output_fol = string(sherlock_fol, sherlock_output_file, input_version,"_",submitdate,"/")
else
    output_fol = string(laptop_fol, laptop_output_file, input_version,"_",submitdate,"/")
end
subsel_data_fol = string(input_fol,inputs[1,:timeseriesID],"/")

# setup unique ID for period identification in output
if multiTF
    pbits = split(periodID,r"_|\.|-")
    periodnum = pbits[2]
    periodfirst = pbits[3]
    periodlast = pbits[4]
    periodsave = string("p",periodnum,"_",periodfirst,"_",periodlast)
end

# --------------------------------------------------------------------------------------
### TIME PERIOD DATA
# -------------------------------------------
if multiTF
    hours = CSV.read(string(subsel_data_fol,periodID),datarow=1,types=[Int])[1]
else
    hours = CSV.read(string(subsel_data_fol,"periods_",timeseriesID,".csv"),datarow=1,types=[Int])[1]
end
n_days = convert(Int32,ceil(length(hours)/24))
n_periods = length(hours)

# convert first and notfirst periods to indices
t_firsts = 1
t_notfirst = collect(2:length(hours))

dem2 = CSV.read(string(default_data_fol, inputs[1,:demandFile]),datarow=2,missingstring="NA")
# subselect for just the rows corresponding to 'hours'
dem = dem2[hours,2]

# --------------------------------------------------------------------------------------
# STOCHASTIC VARIABLE DATA
# -------------------------------------------

## Net Demand uncertainty ##

if !isfile(string(subsel_data_fol,"demandScenarios_vdem","_",inputs[1,:stochID],"_",periodsave,".csv"))
    ndprobs = CSV.read(string(default_data_fol , "dist_input_",inputs[1,:stochID],"_nd.csv"))
    ndv_in = convert(Array,ndprobs[1,:]) # converts the first row of probs to an Array
    ndpro_in = rationalize.(convert(Array,ndprobs[2,:])) #to avoid rounding issues later

    demandScenarios = make_scenarios(n_periods, ndv_in, ndpro_in, int_length; randsel = randScenarioSel, nrand = nrandp)
    vdem = demandScenarios[1]
    pro = demandScenarios[2]
    writecsvmulti(DataFrame(vdem),subsel_data_fol,string("demandScenarios_vdem","_",inputs[1,:stochID]),multiTF,periodsave)
    writecsvmulti(DataFrame(pro),subsel_data_fol,string("demandScenarios_prob","_",inputs[1,:stochID]),multiTF,periodsave)
else
    vdem = convert(Array,CSV.read(string(subsel_data_fol,"demandScenarios_vdem","_",inputs[1,:stochID],"_",periodsave,".csv")))
    pro = convert(Array,CSV.read(string(subsel_data_fol,"demandScenarios_prob","_",inputs[1,:stochID],"_",periodsave,".csv")))
end

# testing:
# if sum(pro) != 1
    # error("sum of probabilities is not one, it is ", sum(pro))
# end
println("sum of probabilities is ", sum(pro))
n_omega = length(pro) #redefine for new number of scenarios

# Set up demand realizations matrix dem_real[t,o]
dem_real = dem .* vdem

# --------------------------------------------------------------------------------------
# GENERATOR DATA
# -------------------------------------------

genset = CSV.read(string(default_data_fol,inputs[1,:genFile]), missingstring ="NA")

# Tips for manipulating genset data:
# genset[Symbol("Plant Name")] # this is how to access by column name if there are spaces
# names(genset) # this is how to get the column names
# anscombe[:,[:X3, :Y1]]  #how to grab several columns by colname

dr_ind = findin(genset[:Fuel],dr_gens)
slow_ind = findin(genset[:Speed],["SLOW"])
fast_ind = findin(genset[:Speed],["FAST"])

# adjust slow/fast defn based on DRtype
if DRtype == 1
    fast_ind = vcat(fast_ind,dr_ind)
elseif (DRtype == 2) | (DRtype == 3)
    slow_ind = vcat(slow_ind,dr_ind)
elseif DRtype == 0
    println("No DR in this run")
else
    error("DRtype must be 0, 1, 2 or 3")
end

#generator min and max
pmin = genset[:PMin]
pmax = genset[:Capacity]
startup = genset[:StartCost]
varcost = genset[:VCost]
rampmax = genset[:ramprate]

# for manual override of DR variable cost
if dr_override
    varcost[dr_ind] = dr_varcost
end

# Generator PARAMS ###
n_gsl= length(slow_ind)# number of slow generators
n_g =nrow(genset)# number of generators
n_gf = length(fast_ind)
n_gdr = length(dr_ind) #number of DR generators
n_t = n_periods # number of timesteps

## GENERATOR AVAILABILITY ##
pf = repeat([1.0], inner = [n_g, n_t])

### Wind and solar
# load wind, solar info
solar_avail = CSV.read(string(default_data_fol,"solar_availability_factors_2016.csv"))
wind_avail = CSV.read(string(default_data_fol,"wind_availability_factors_2016.csv"))
# remove last col of each
solar_avail = solar_avail[:,1:(ncol(solar_avail)-1)]
wind_avail = wind_avail[:,1:(ncol(wind_avail)-1)]
# loop through all colnames, use findin(genset[:plantUnique],XX) to get row
# sub in new info
for i in 1:length(names(solar_avail))
    col = names(solar_avail)[i]
    ind = findin(genset[:plantUnique],[convert(String, col)])
    pf[ind,:] = solar_avail[hours,i]
end

for i in 1:length(names(wind_avail))
    col = names(wind_avail)[i]
    ind = findin(genset[:plantUnique],[convert(String, col)])
    pf[ind,:] = wind_avail[hours,i]
end

### Demand Response
if parse(inputs[1,:availID])!=0 &  DRtype != 0
    dr_avail = CSV.read(string(default_data_fol,inputs[1,:availID],".csv"))
    # remove first col
    dr_avail = dr_avail[:,2:ncol(dr_avail)]
    # sub in new info
    for i in 1:length(names(dr_avail))
        col = names(dr_avail)[i]
        ind = findin(genset[:plantUnique],[convert(String, col)])
        if length(ind) > 0
            pf[ind,:] = dr_avail[hours,i]
        else
            error("No corresponding DR unit found for ", names(dr_avail)[i])
        end
    end
end

# --------------------------------------------------------------------------------------
# MODEL SETUP
# -------------------------------------------
### SETS ###
TIME = 1:n_t
SCENARIOS = 1:n_omega
GENERATORS = 1:n_g #all generators
GEN_NODR = findin(genset[:Fuel],notdr_gens)
GF = fast_ind
GSL = slow_ind #slow generators
GDR = dr_ind #DR generators
# index for where the DR is in the slow generators
GDR_SL_ind = findin(GSL,GDR)

# -------------------------------------------
### MODEL ###
# m = Model(solver = ClpSolver())
m = Model(solver=GurobiSolver(Method=MethodParam,#1,
                              MIPFocus=MIPFocusParam,#3,
                              MIPGap=MIPGapParam,#0.0003,
                              NodefileStart = NodefileStartParam,#0.05,
                              NodefileDir = "/scratch/users/pjlevi/gurobi_solving_outputs/",
                              Seed = convert(Int64,abs(floor(rand(Float64)*2000000000)))))

if no_vars
    error("just testing model so we are stopping here")
else
    println("begin building model")
end

# -----------------------------------------------------------------------------------------------------
### VARIABLES ###
# @variable(m, 0 <= x <= 2 )
@variable(m, z[1:n_gsl,1:n_t]) # slow generator startup
if trueBinaryStartup
    @variable(m, w[1:n_gsl,1:n_t], Bin) # slow generator commitment TRUE BINARY
else
    @variable(m, 0 <= w[1:n_gsl,1:n_t] <= 1) # slow generator commitment RELAXED BINARY
end

#real-time commitment, startup
@variable(m, v[1:n_g,1:n_t,1:n_omega]) # generator startup
if trueBinaryStartup
    @variable(m, u[1:n_g,1:n_t,1:n_omega], Bin) # generator commitment TRUE BINARY
else
    @variable(m, 0 <= u[1:n_g,1:n_t,1:n_omega] <= 1) # generator commitment RELAXED BINARY
end

#production variables
@variable(m, p[1:n_g,1:n_t,1:n_omega] >= 0) #generator production

# if DR requires advance production schedule
if DRtype == 3
    @variable(m, p_dr[1:n_gdr, 1:n_t] >= 0) # DR day-ahead production commitment
end

# track number of startups to determine costs
@variable(m, start_num[1:n_g, 1:n_t, 1:n_omega] >= 0)

# -----------------------------------------------------------------------------------------------------
### CONSTRAINTS ###

#SUPPLY-DEMAND
@constraint(m,supplydemand[t=1:n_t, o=1:n_omega],
    sum( p[g,t,o] for g=1:n_g) >= dem_real[t,o])

#GENMIN
@constraint(m, mingen[g= 1:n_g, t= 1:n_t, o=1:n_omega ],
    p[g,t,o] >= pmin[g] * u[g,t,o])
#GENMAX - only for GEN_NODR if dr is random
@constraint(m,[g = 1:n_g, t = 1:n_t, o=1:n_omega],
    p[g,t,o] <= pmax[g] * u[g,t,o] * pf[g,t] )

if DRtype == 3
#GENMAXDR
    @constraint(m,[g = 1:n_gdr, t = 1:n_t],
        p_dr[g,t] <= pmax[GDR[g]] * w[GDR_SL_ind[g],t] * pf[GDR_SL_ind[g],t]) #needed for p_da
    @constraint(m,type3dr[g=1:n_gdr,t = TIME, o = SCENARIOS],
         p[GDR[g],t,o] == p_dr[g,t])
end

#START_S
@constraint(m,[g = 1:n_gsl, t=t_notfirst],
    z[g,t] == w[g,t] - w[g,t-1])
#START_F
@constraint(m,[g=GENERATORS ,t=t_notfirst, o = SCENARIOS],
    v[g,t,o] == u[g,t,o] - u[g,t-1,o])
#INIT_S
@constraint(m,[g=1:n_gsl,t=t_firsts],
    z[g,t] == w[g,t])
#INIT_F
@constraint(m,[g=GENERATORS,t=t_firsts, o = SCENARIOS],
    v[g,t,o] == u[g,t,o])
#NAN_ST
@constraint(m,[g = 1:n_gsl, t = TIME, o = SCENARIOS],
    v[GSL[g],t,o] == z[g,t])
#NAN_CM
@constraint(m,[g = 1:n_gsl, t = TIME, o = SCENARIOS],
    u[GSL[g],t,o] == w[g,t])

#STARTUP COUNT
@constraint(m, [g=GENERATORS, t = TIME, o = SCENARIOS],
    start_num[g,t,o] >= v[g,t,o])

#RAMP RATE -
if ramplims !=0
    @constraint(m,ramplim[g=GENERATORS,t=t_notfirst, o = SCENARIOS],
        p[g,t,o] - p[g,t-1,o] <= (rampmax[g] * pmax[g]))
end

# ------ DR USAGE LIMITS

#STARTLIM
# number of startups per period
if startlim !=0
    # # per day
    # @constraint(m,[g = 1:n_gdr,d = 1:n_days, o = SCENARIOS],
    #     sum(v[GDR[g],t,o] for t = (24*(d-1)+1):(24*d)) <= startlim)
    # per period
    @constraint(m,limstart[g = 1:n_gdr, o = SCENARIOS],
        sum(start_num[GDR[g],t,o] for t = TIME) <= startlim)
end

#HOURLIM
#number of hours used per period
if hourlim != 0
    # # per day
    # @constraint(m,[g = 1:n_gdr,d = 1:n_days, o = SCENARIOS],
    #     sum(u[GDR[g],t,o] for t = (24*(d-1)+1):(24*d)) <= hourlim)
    # per period
    @constraint(m,limhrs[g = 1:n_gdr, o = SCENARIOS],
        sum(u[GDR[g],t,o] for t = TIME) <= hourlim)
end

#ENERGYLIM
# amount of energy used per period
if energylim != 0
    # per day
    # @constraint(m,[g = 1:n_gdr,d = 1:n_days ,o=SCENARIOS],
    #     sum(p[GDR[g],t,o] for t = (24*(d-1)+1):(24*d)) <= energylim)
    # per period
    @constraint(m,limenergy[g = 1:n_gdr, o=SCENARIOS],
        sum(p[GDR[g],t,o] for t = TIME) <= energylim)
end

#DURATIONLIM
# number of hours in a row that DR can be committed
if durationlim != 0
    # make a list of lists of 24-h timestep groupings
    l = Int(durationlim + 1)
    s = 1
    WINDOWS = [TIME[i:i+l-1] for i in range(1,s,Int(length(TIME)-l+1))]
    for i in range(1,s,Int(length(TIME)-l+1))
        print(TIME[i:i+l-1])
    end
    for W in WINDOWS
        @constraint(m,[g = 1:n_gdr, o = SCENARIOS],
            sum(u[GDR[g],t,o] for t = W) <= durationlim)
    end
end

### OBJECTIVE ###
@objective(m, Min, sum(pro[o] *
    sum(start_num[g,t,o]*startup[g] + p[g,t,o]*varcost[g] for g = GENERATORS, t = TIME)
    for o = SCENARIOS))

# Check model
# print(m)
if debug
    error("just testing model so we are stopping here")
end

# -------------------------------------------
# SOLVE MODEL
# -------------------------------------------
@printf("\nSolving:\n")
status = solve(m)
@printf("Status: %s\n", status)

if !isdir(output_fol)
    mkdir(output_fol)
end

## save a copy of inputs & demand scenarios to the output file
writecsvmulti(read_inputs,output_fol,"inputfile",multiTF,periodsave)

# --------------------------------------------------------------------------------------
# SAVE OUTPUT
# -------------------------------------------

# production of DR
y = getvalue(p[GDR,:,:])

if DRtype != 0
    y_out = convert3dto2d(y,1, 3,  2,
        vcat([String("o$i") for i in 1:n_omega],"DR_IND","t"),
         genset[dr_ind,:plantUnique])
    writecsvmulti(y_out,output_fol,"DR_production",multiTF,periodsave)
end

# production of slow generators
zs = getvalue(p[GSL,:,:])

y_out = convert3dto2d(zs,1, 3, 2,
    vcat([String("o$i") for i in 1:n_omega],"GEN_IND","t"),
     genset[slow_ind,:plantUnique])
writecsvmulti(y_out,output_fol,"slow_production",multiTF,periodsave)

# production of fast generators
zf = getvalue(p[GF,:,:])
y_out = convert3dto2d(zf,1, 3, 2,
    vcat([String("o$i") for i in 1:n_omega],"GEN_IND","t"),
     genset[fast_ind,:plantUnique])
writecsvmulti(y_out,output_fol,"fast_production",multiTF,periodsave)

# commitment of slow generators
w_out = getvalue(w)
wdf = DataFrame(transpose(w_out))
names!(wdf,[Symbol("$input") for input in genset[slow_ind,:plantUnique]])
writecsvmulti(wdf,output_fol,"slow_commitment",multiTF,periodsave)

# commitment of all generators
u_out = getvalue(u)
y_out = convert3dto2d(u_out,1, 3, 2,
    vcat([String("o$i") for i in 1:n_omega],"GEN_IND","t"),
     genset[:plantUnique])
writecsvmulti(y_out,output_fol,"u_commitment",multiTF,periodsave)

# startup of all generators
v_out = getvalue(v)
y_out = convert3dto2d(v_out,1, 3, 2,
    vcat([String("o$i") for i in 1:n_omega],"GEN_IND","t"),
     genset[:plantUnique])
writecsvmulti(y_out,output_fol,"v_startup",multiTF,periodsave)

# costs
totcost = getvalue(sum(pro[o] *
    sum(start_num[g,t,o]*startup[g] + p[g,t,o]*varcost[g] for g = GENERATORS, t = TIME)
    for o = SCENARIOS))
totstartupcost = getvalue(sum(pro[o] *
    sum(start_num[g,t,o]*startup[g] for g = GENERATORS, t = TIME)
    for o = SCENARIOS))
totvarcost = getvalue(sum(pro[o] *
    sum(p[g,t,o]*varcost[g] for g = GENERATORS, t = TIME)
    for o = SCENARIOS))

output_summary = DataFrame(TotalCost = totcost, TotStartupCst = totstartupcost,
                            TotVarCst = totvarcost,
                            FracStartupCost = totstartupcost/totcost,
                            FracVarCost = totvarcost/totcost,
                            timeseries = inputs[1,:timeseriesID],
                            stoch_scenario = inputs[1,:stochID])

writecsvmulti(output_summary,output_fol,"summary_stats",multiTF,periodsave)

# Get dual variables and save
# only works if the problem is linear - cannot have binary variables
if !trueBinaryStartup
    #supplydemand #2D
    sd_shadow = DataFrame(getdual(supplydemand))
    writecsvmulti(sd_shadow,output_fol,"supplydemand_shadow",multiTF,periodsave)

    ### mingen ###
    mingen_shadow = getdual(mingen)
    mingen_sf = convert3dto2d(mingen_shadow,1, 3,  2,
        vcat([String("o$i") for i in 1:n_omega],"GEN_IND","t"),
         genset[:,:plantUnique])
    # only save rows where there is a nonzero shadow price
    shadowsum = sum(convert(Array{Float64},mingen_sf[:,1:n_omega]),2)
    saveind = find(shadowsum .!= 0)
    writecsvmulti(mingen_sf[saveind,:],output_fol,"mingen_shadow",multiTF,periodsave)

    ### type3dr ###
    ## TODO: This does not work.
    if DRtype == 3
        type3dr_shadow = getdual(type3dr)
        type3dr_sf = convert3dto2d(type3dr_shadow,1, 3,  2,
            vcat([String("o$i") for i in 1:n_omega],"DR_IND","t"),
             genset[dr_ind,:plantUnique])
        # only save rows where there is a nonzero shadow price
        shadowsum = sum(convert(Array{Float64},type3dr_sf[:,1:n_omega]),2)
        saveind = find(shadowsum .!= 0)
        writecsvmulti(type3dr_sf[saveind,:],output_fol,"type3dr_shadow",multiTF,periodsave)
    end
end

max_stars_repo_name: pajalevi/Stoch-UC-model-for-DR

max_issues_repo_name: pajalevi/Stoch-UC-model-for-DR

max_forks_repo_name: pajalevi/Stoch-UC-model-for-DR


text:
# NMF for dense matrices

using NMF

function run(algname)

    # prepare data
    p = 8
    k = 5
    n = 100

    Wg = abs(randn(p, k))
    Hg = abs(randn(k, n))
    X = Wg * Hg + 0.1 * randn(p, n)

    # run NNMF
    println("Algorithm: $(algname)")
    println("---------------------------------")

    r = nnmf(X, k;
             init=:nndsvdar,
             alg=symbol(algname),
             maxiter=30,
             verbose=true)

    # display results
    println("numiters  = $(r.niters)")
    println("converged = $(r.converged)")
    @printf("objvalue  = %.6e\n", r.objvalue)
    println("W matrix = ")
    NMF.printf_mat(r.W)

    println()
end

function print_help()
    println("Usage:")
    println()
    println("  julia densenmf.jl <alg>")
    println()
    println("  <alg> is the name of the chosen algorithm, which can be ")
    println()
    println("    multmse:   Multiplicative update (minimize MSE)")
    println("    multdiv:   Multiplicative update (minimize divergence)")
    println("    projals:   Projected ALS")
    println("    alspgrad:  ALS Projected Gradient Descent")
    println()
end

function main(args)
    if isempty(args)
        print_help()
    elseif length(args) == 1
        a = lowercase(args[1])
        if a == "-h" || a == "--help"
            print_help()
        else
            run(a)
        end
    else
        @warn("Invalid command line arguments.")
        print_help()
    end
end

main(ARGS)


text:
push!(LOAD_PATH, "../src/")

using Documenter
using PointProcesses

DocMeta.setdocmeta!(
    PointProcesses,
    :DocTestSetup,
    :(using PointProcesses);
    recursive = true,
)

makedocs(;
    modules = [PointProcesses],
    authors = "Guillaume Dalle",
    repo = "https://github.com/gdalle/PointProcesses.jl/blob/{commit}{path}#{line}",
    sitename = "PointProcesses.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://gdalle.github.io/PointProcesses.jl",
        assets = String[],
    ),
    pages = [
        "Home" => "index.md",
        "Tutorial" => [
            "Event history" => "tuto/history.md",
            "Point processes" => "tuto/point_processes.md",
            "Built-in models" => "tuto/models.md",
            "Markov chains" => "tuto/markov.md",
            "Hidden Markov models" => "tuto/hmm.md",
            "Utilities" => "tuto/utils.md",
        ],
        "API reference" => "api.md",
        "Roadmap" => "roadmap.md",
    ],
)

deploydocs(; repo = "github.com/gdalle/PointProcesses.jl", devbranch = "master")


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

function restore!(model::MOI.ModelLike, ::AbstractIntegerAlgorithm)
    return nothing
end

function check_optimality(lshaped::AbstractLShaped, ::AbstractIntegerAlgorithm)
    return all(1:num_subproblems(lshaped)) do idx
        if has_metadata(lshaped.execution.metadata, idx, :integral_solution)
            return get_metadata(lshaped.execution.metadata, idx, :integral_solution)
        else
            error("Worker $idx has not reported if latest solution satisfies integral restrictions.")
        end
    end
end

function is_approx_integer(val::AbstractFloat)
    return ceil(val) - val <= sqrt(eps()) || val - floor(val) <= sqrt(eps())
end

function gather_integer_variables(subproblem::SubProblem)
    model = subproblem.model
    N = num_stages(model.ext[:decisions])
    all_known = mapreduce(vcat, 1:N-1) do s
        index.(all_known_decision_variables(model, s))
    end
    all_decisions = index.(all_decision_variables(model, N))
    indices = Vector{MOI.VariableIndex}()
    for var in all_variables(model)
        vi = index(var)
        if vi in all_known
            # Known decision, skip
            continue
        end
        if vi in all_decisions
            # Decision variable
            dref = DecisionRef(model, vi)
            if is_integer(dref) || is_binary(dref)
                push!(indices, vi)
            else
                continue
            end
        else
            if is_integer(var) || is_binary(var)
                push!(indices, vi)
            else
                continue
            end
        end
    end
    return indices
end

function check_integrality_restrictions(subproblem::SubProblem, fractional = Vector{MOI.VariableIndex}())
    satisfies_integer = map(integer_variables(subproblem.integer_algorithm)) do vi
        val = MOI.get(subproblem.optimizer, MOI.VariablePrimal(), vi)
        if is_approx_integer(MOI.get(subproblem.optimizer, MOI.VariablePrimal(), vi))
            return true
        else
            push!(fractional, vi)
            return false
        end
    end
    return all(satisfies_integer)
end

mutable struct SecondStageLP{T_ <: AbstractFloat}
    T::SparseMatrixCSC{T_,Int}
    W::SparseMatrixCSC{T_,Int}
    h::Vector{T_}
    lb::Vector{T_}
    ub::Vector{T_}
end

function SecondStageLP(subproblem::SubProblem{T_}; standard_form = false) where T_ <: AbstractFloat
    model = subproblem.model
    nx = StochasticPrograms.num_known_decisions(model)
    ny = num_variables(model) - nx
    # Prepare data
    Tᵢ = Vector{Int}()
    Tⱼ = Vector{Int}()
    Tᵥ = Vector{T_}()
    Wᵢ = Vector{Int}()
    Wⱼ = Vector{Int}()
    Wᵥ = Vector{T_}()
    lhs = Vector{T_}()
    rhs = Vector{T_}()
    lb = fill(-Inf, ny)
    ub = fill(Inf, ny)
    i = 0
    slack = 0
    # Loop over constraints
    for (F,S) in list_of_constraint_types(model)
        for cref in all_constraints(model, F, S)
            if S <: StochasticPrograms.SingleDecisionSet ||
                S <: StochasticPrograms.MultipleDecisionSet
                continue
            elseif F <: VariableRef || F <: DecisionRef
                f = MOI.get(model, MOI.ConstraintFunction(), cref)
                if f isa MOI.SingleVariable
                    idx = f.variable.value - nx
                elseif f isa SingleDecision
                    idx = f.decision.value - nx
                end
                set = MOI.get(model, MOI.ConstraintSet(), cref)
                if set isa MOI.LessThan
                    ub[idx] = set.upper
                    if standard_form
                        i += 1
                        slack += 1
                        push!(Wᵢ, i)
                        push!(Wⱼ, idx)
                        push!(Wᵥ, 1.0)
                        push!(Wᵢ, i)
                        push!(Wⱼ, ny + slack)
                        push!(Wᵥ, 1.0)
                        push!(lhs, -Inf)
                        push!(rhs, set.upper)
                    end
                elseif set isa MOI.GreaterThan
                    lb[idx] = set.lower
                elseif set isa MOI.EqualTo
                    lb[idx] = set.value
                    ub[idx] = set.value
                end
            elseif F <: AffExpr
                i += 1
                f = MOI.get(model, MOI.ConstraintFunction(), cref)
                for term in f.terms
                    push!(Wᵢ, i)
                    push!(Wⱼ, term.variable_index.value - nx)
                    push!(Wᵥ, term.coefficient)
                end
                set = MOI.get(model, MOI.ConstraintSet(), cref)
                if set isa MOI.LessThan
                    push!(lhs, -Inf)
                    push!(rhs, set.upper)
                    if standard_form
                        slack += 1
                        push!(Wᵢ, i)
                        push!(Wⱼ, ny + slack)
                        push!(Wᵥ, 1.0)
                    end
                elseif set isa MOI.GreaterThan
                    push!(lhs, set.lower)
                    push!(rhs, Inf)
                    if standard_form
                        slack += 1
                        push!(Wᵢ, i)
                        push!(Wⱼ, ny + slack)
                        push!(Wᵥ, -1.0)
                    end
                elseif set isa MOI.EqualTo
                    push!(lhs, set.value)
                    push!(rhs, set.value)
                    if standard_form
                        slack += 1
                        push!(Wᵢ, i)
                        push!(Wⱼ, ny + slack)
                        push!(Wᵥ, 1.0)
                    end
                end
            elseif F <: Vector{<:AffExpr}
                f = MOI.get(model, MOI.ConstraintFunction(), cref)
                set = MOI.get(model, MOI.ConstraintSet(), cref)
                for (j,fⱼ) in enumerate(MOIU.eachscalar(f))
                    i += 1
                    for term in fⱼ.terms
                        push!(Wᵢ, i)
                        push!(Wⱼ, term.variable_index.value - nx)
                        push!(Wᵥ, term.coefficient)
                    end
                    if set isa MOI.Nonpositives
                        push!(lhs, -Inf)
                        push!(rhs, -f.constants[j])
                        if standard_form
                            slack += 1
                            push!(Wᵢ, i)
                            push!(Wⱼ, ny + slack)
                            push!(Wᵥ, 1.0)
                        end
                    elseif set isa MOI.Nonnegatives
                        push!(lhs, -f.constants[j])
                        push!(rhs, Inf)
                        if standard_form
                            slack += 1
                            push!(Wᵢ, i)
                            push!(Wⱼ, ny + slack)
                            push!(Wᵥ, -1.0)
                        end
                    elseif set isa MOI.Zeros
                        push!(lhs, -f.constants[j])
                        push!(rhs, -f.constants[j])
                        if standard_form
                            slack += 1
                            push!(Wᵢ, i)
                            push!(Wⱼ, ny + slack)
                            push!(Wᵥ, 1.0)
                        end
                    end
                end
            elseif F <: DecisionAffExpr
                i += 1
                aff = jump_function(model, MOI.get(model, MOI.ConstraintFunction(), cref))
                for (coeff, var) in linear_terms(aff.variables)
                    push!(Wᵢ, i)
                    push!(Wⱼ, index(var).value - nx)
                    push!(Wᵥ, coeff)
                end
                for (coeff, dvar) in linear_terms(aff.decisions)
                    if state(dvar) == Known
                        push!(Tᵢ, i)
                        push!(Tⱼ, index(dvar).value)
                        push!(Tᵥ, coeff)
                    else
                        push!(Wᵢ, i)
                        push!(Wⱼ, index(dvar).value - nx)
                        push!(Wᵥ, coeff)
                    end
                end
                set = MOI.get(model, MOI.ConstraintSet(), cref)
                if set isa MOI.LessThan
                    push!(lhs, -Inf)
                    push!(rhs, set.upper)
                    if standard_form
                        slack += 1
                        push!(Wᵢ, i)
                        push!(Wⱼ, ny + slack)
                        push!(Wᵥ, 1.0)
                    end
                elseif set isa MOI.GreaterThan
                    push!(lhs, set.lower)
                    push!(rhs, Inf)
                    if standard_form
                        slack += 1
                        push!(Wᵢ, i)
                        push!(Wⱼ, ny + slack)
                        push!(Wᵥ, -1.0)
                    end
                elseif set isa MOI.EqualTo
                    push!(lhs, set.value)
                    push!(rhs, set.value)
                    if standard_form
                        slack += 1
                        push!(Wᵢ, i)
                        push!(Wⱼ, ny + slack)
                        push!(Wᵥ, 1.0)
                    end
                end
            elseif F <: Vector{<:DecisionAffExpr}
                f = MOI.get(model, MOI.ConstraintFunction(), cref)
                set = MOI.get(model, MOI.ConstraintSet(), cref)
                for (j,fⱼ) in enumerate(MOIU.eachscalar(f))
                    i += 1
                    aff = jump_function(model, fⱼ)
                    for (coeff, var) in linear_terms(aff.variables)
                        push!(Wᵢ, i)
                        push!(Wⱼ, index(var).value - nx)
                        push!(Wᵥ, coeff)
                    end
                    for (coeff, dvar) in linear_terms(aff.decisions)
                        if state(dvar) == Known
                            push!(Tᵢ, i)
                            push!(Tⱼ, index(dvar).value)
                            push!(Tᵥ, coeff)
                        else
                            push!(Wᵢ, i)
                            push!(Wⱼ, index(dvar).value - nx)
                            push!(Wᵥ, coeff)
                        end
                    end
                    if set isa MOI.Nonpositives
                        push!(lhs, -Inf)
                        push!(rhs, -f.variable_part.constants[j])
                        if standard_form
                            slack += 1
                            push!(Wᵢ, i)
                            push!(Wⱼ, ny + slack)
                            push!(Wᵥ, 1.0)
                        end
                    elseif set isa MOI.Nonnegatives
                        push!(lhs, -f.variable_part.constants[j])
                        push!(rhs, Inf)
                        if standard_form
                            slack += 1
                            push!(Wᵢ, i)
                            push!(Wⱼ, ny + slack)
                            push!(Wᵥ, -1.0)
                        end
                    elseif set isa MOI.Zeros
                        push!(lhs, -f.variable_part.constants[j])
                        push!(rhs, -f.variable_part.constants[j])
                        if standard_form
                            slack += 1
                            push!(Wᵢ, i)
                            push!(Wⱼ, ny + slack)
                            push!(Wᵥ, 1.0)
                        end
                    end
                end
            else
                error("Cannot extract linear representation from $F-$S constraint.")
            end
        end
    end
    T̃ = sparse(Tᵢ, Tⱼ, Tᵥ, i, nx)
    W̃ = sparse(Wᵢ, Wⱼ, Wᵥ, i, ny + slack)
    if standard_form
        h = lhs
        i = .!isfinite.(h)
        h[i] .= rhs[i]
        return SecondStageLP(T̃, W̃, h, lb, ub)
    else
        T, h = StochasticPrograms.canonical(T̃, lhs, rhs)
        W, h = StochasticPrograms.canonical(W̃, lhs, rhs)
        return SecondStageLP(-T, -W, -h, lb, ub)
    end
end

function optimal_basis(subproblem::SubProblem)
    model = subproblem.model
    nx = StochasticPrograms.num_known_decisions(model)
    ny = num_variables(model) - nx
    variable_basis = Int[]
    slack_basis = Int[]
    slack = 0
    # Loop over constraints
    for (F,S) in list_of_constraint_types(model)
        for cref in all_constraints(model, F, S)
            if S <: StochasticPrograms.SingleDecisionSet ||
                S <: StochasticPrograms.MultipleDecisionSet
                continue
            elseif F <: VariableRef || F <: DecisionRef
                f = MOI.get(model, MOI.ConstraintFunction(), cref)
                set = MOI.get(model, MOI.ConstraintSet(), cref)
                if f isa MOI.SingleVariable
                    idx = f.variable.value - nx
                elseif f isa SingleDecision
                    idx = f.decision.value - nx
                end
                if set isa MOI.LessThan || set isa MOI.EqualTo
                    slack += 1
                    if MOI.get(model, MOI.ConstraintBasisStatus(), cref) == MOI.BASIC
                        push!(slack_basis, ny + slack)
                    end
                elseif set isa MOI.GreaterThan
                    if MOI.get(model, MOI.ConstraintBasisStatus(), cref) == MOI.BASIC
                        push!(variable_basis, idx)
                    end
                end
            elseif F <: AffExpr || F <: DecisionAffExpr
                slack += 1
                if MOI.get(model, MOI.ConstraintBasisStatus(), cref) == MOI.BASIC
                    push!(slack_basis, ny + slack)
                end
            else
                error("Cannot extract basis from $F-$S constraint.")
            end
        end
    end
    return vcat(variable_basis, slack_basis)
end

function add_row!(lp::SecondStageLP, T̄::AbstractVector, W̄::AbstractVector, h̄::AbstractFloat; standard_form = false)
    # Get current size
    m = size(lp.W, 1)
    nx = size(lp.T, 2)
    ny = size(lp.W, 2)
    if standard_form
        # Prepare to add new slack variable
        ny += 1
    end
    # Get current matrix structure
    Tᵢ, Tⱼ, Tᵥ = findnz(lp.T)
    Wᵢ, Wⱼ, Wᵥ = findnz(lp.W)
    # Add T row
    for (j, coeff) in enumerate(T̄)
        if abs(coeff) >= sqrt(eps())
            push!(Tᵢ, m + 1)
            push!(Tⱼ, j)
            push!(Tᵥ, coeff)
        end
    end
    lp.T = sparse(Tᵢ, Tⱼ, Tᵥ, m + 1, nx)
    # Add W row
    for (j, coeff) in enumerate(W̄)
        if abs(coeff) >= sqrt(eps())
            push!(Wᵢ, m + 1)
            push!(Wⱼ, j)
            push!(Wᵥ, coeff)
        end
    end
    if standard_form
        # Add slack term
        push!(Wᵢ, m + 1)
        push!(Wⱼ, ny)
        push!(Wᵥ, -1.0)
    end
    lp.W = sparse(Wᵢ, Wⱼ, Wᵥ, m + 1, ny)
    # Add new rhs
    push!(lp.h, h̄)
    return nothing
end


text:
module CMUDict

export cmudict

using ..Data: deps

const version = "0.7b"

function load()
  isdir(deps("cmudict")) && return
  mkpath(deps("cmudict"))
  for x in ["", ".phones", ".symbols"]
    download("http://svn.code.sf.net/p/cmusphinx/code/trunk/cmudict/cmudict-$version$x",
             deps("cmudict", "cmudict$x"))
  end
end

function phones()
  load()
  Symbol.(first.(split.(split(readstring(deps("cmudict", "cmudict.phones")),
                        "\n", keep = false), "\t")))
end

function symbols()
  load()
  Symbol.(split(readstring(deps("CMUDict", "cmudict.symbols")),
                "\n", keep = false))
end

function rawdict()
  load()
  Dict(String(xs[1]) => Symbol.(xs[2:end]) for xs in
       filter(!isempty, split.(split(readstring(deps("CMUDict", "cmudict")), "\n"))))
end

validword(s) = ismatch(r"^[\w-\.]+$", s)

cmudict() = filter((s, ps) -> validword(s), rawdict())

alphabet() = ['A':'Z'..., '0':'9'..., '_', '-', '.']

end


text:
Base.promote_rule{T1<:AbstractGray,T2<:AbstractGray}(::Type{T1}, ::Type{T2}) = Gray{promote_type(eltype(T1), eltype(T2))}
Base.promote_rule{T1<:AbstractRGB,T2<:AbstractRGB}(::Type{T1}, ::Type{T2}) = RGB{promote_type(eltype(T1), eltype(T2))}

# no-op and element-type conversions, plus conversion to and from transparency
# Colorimetry conversions are in Colors.jl
convert{C<:Colorant}(::Type{C}, c::C) = c
convert{C<:Colorant}(::Type{C}, c) = cconvert(ccolor(C, typeof(c)), c)
cconvert{C}(::Type{C}, c::C) = c
cconvert{C}(::Type{C}, c)    = _convert(C, base_color_type(C), base_color_type(c), c)
convert{C<:TransparentColor}(::Type{C}, c::Color, alpha) = cconvert(ccolor(C, typeof(c)), c, alpha)
cconvert{C<:Color,T,N}(::Type{AlphaColor{C,T,N}}, c::C, alpha) = alphacolor(C)(c, alpha)
cconvert{C<:Color,T,N}(::Type{ColorAlpha{C,T,N}}, c::C, alpha) = coloralpha(C)(c, alpha)
cconvert{C<:TransparentColor}(::Type{C}, c::Color, alpha) =_convert(C, base_color_type(C), base_color_type(c), c, alpha)

# Fallback definitions that print nice error messages
_convert{C}(::Type{C}, ::Any, ::Any, c) = error("No conversion of ", c, " to ", C, " has been defined")
_convert{C}(::Type{C}, C1::Any, C2::Any, c, alpha) = error("No conversion of (", c, ",alpha=$alpha) to ", C, " with consistency-types $C1 and $C2 has been defined")

# Any AbstractRGB types can be interconverted
# (the first 2 are just for ambiguity resolution)
# Note: on julia 0.3 these have to be before the block below, or you
# get a spurious ambiguity warning.
_convert{Cout<:AbstractRGB,C1<:AbstractRGB}(::Type{Cout}, ::Type{C1}, ::Type{C1}, c) = Cout(red(c), green(c), blue(c))
_convert{A<:TransparentRGB,C1<:AbstractRGB}(::Type{A}, ::Type{C1}, ::Type{C1}, c) = A(red(c), green(c), blue(c), alpha(c))
_convert{Cout<:AbstractRGB,C1<:AbstractRGB,C2<:AbstractRGB}(::Type{Cout}, ::Type{C1}, ::Type{C2}, c) = Cout(red(c), green(c), blue(c))
_convert{A<:TransparentRGB,C1<:AbstractRGB,C2<:AbstractRGB}(::Type{A}, ::Type{C1}, ::Type{C2}, c) = A(red(c), green(c), blue(c), alpha(c))

# Implementations for when the base color type is not changing
# These might trip/add transparency, however
_convert{Cout<:Color3,Ccmp<:Color3}(::Type{Cout}, ::Type{Ccmp}, ::Type{Ccmp}, c) = Cout(comp1(c), comp2(c), comp3(c))
_convert{A<:Transparent3,Ccmp<:Color3}(::Type{A}, ::Type{Ccmp}, ::Type{Ccmp}, c) = A(comp1(c), comp2(c), comp3(c), alpha(c))
_convert{Cout<:AbstractGray,Ccmp<:AbstractGray}(::Type{Cout}, ::Type{Ccmp}, ::Type{Ccmp}, c) = Cout(gray(c))
_convert{A<:TransparentGray,Ccmp<:AbstractGray}(::Type{A}, ::Type{Ccmp}, ::Type{Ccmp}, c) = A(gray(c), alpha(c))

# With user-supplied alpha
_convert{A<:Transparent3,Ccmp<:Color3}(::Type{A}, ::Type{Ccmp}, ::Type{Ccmp}, c, alpha) = A(comp1(c), comp2(c), comp3(c), alpha)

# Grayscale
_convert{Cout<:AbstractGray,C1<:AbstractGray,C2<:AbstractGray}(::Type{Cout}, ::Type{C1}, ::Type{C2}, c) = Cout(gray(c))
_convert{A<:TransparentGray,C1<:AbstractGray,C2<:AbstractGray}(::Type{A}, ::Type{C1}, ::Type{C2}, c) = A(gray(c), alpha(c))
_convert{A<:TransparentGray,C1<:AbstractGray,C2<:AbstractGray}(::Type{A}, ::Type{C1}, ::Type{C2}, c, alpha) = A(gray(c), alpha)

for T in (RGB24, ARGB32, Gray24, AGray32)
    @eval begin
        @deprecate convert(::Type{UInt32}, c::$T) reinterpret(UInt32, c)
        @deprecate getindex(::Type{$T}, x::UInt32) $T[reinterpret($T, x)]
        @deprecate getindex(::Type{$T}, x::UInt32, y::UInt32) $T[reinterpret($T, x), reinterpret($T, y)]
        @deprecate getindex(::Type{$T}, x::UInt32, y::UInt32, z::UInt32) $T[reinterpret($T, x), reinterpret($T, y), reinterpret($T, z)]
        @deprecate getindex(::Type{$T}, vals::UInt32...) $T[map(x->reinterpret($T, x), vals)...]
    end
end

convert(::Type{RGB24},   x::Real) = RGB24(x, x, x)
convert(::Type{ARGB32},  x::Real) = ARGB32(x, x, x, 1)
convert(::Type{Gray24},  x::Real) = Gray24(x)
convert(::Type{AGray32}, x::Real) = AGray32(x, 1)
convert(::Type{AGray32}, x::Real, alpha) = AGray32(x, alpha)

convert{T}(::Type{Gray{T}},  x::Real)    = Gray{T}(x)
convert{T<:Real}(::Type{T},  x::Gray)    = convert(T, x.val)
convert{T<:Real}(::Type{T},  x::Gray24)  = convert(T, gray(x))
convert{T}(::Type{AGray{T}}, x::Real)    = AGray{T}(x)
convert{T}(::Type{GrayA{T}}, x::Real)    = GrayA{T}(x)

# Generate the transparent analog of a color
alphacolor{C<:Color}(c::C) = alphacolor(C)(c)
alphacolor{C<:Color}(c::C,a) = alphacolor(C)(c,a)
alphacolor{C<:TransparentColor}(c::C) = alphacolor(base_color_type(C))(color(c), alpha(c))
alphacolor{C<:TransparentColor}(c::C,a) = alphacolor(base_color_type(C))(color(c), a)
coloralpha{C<:Color}(c::C) = coloralpha(C)(c)
coloralpha{C<:Color}(c::C,a) = coloralpha(C)(c,a)
coloralpha{C<:TransparentColor}(c::C) = coloralpha(base_color_type(C))(color(c), alpha(c))
coloralpha{C<:TransparentColor}(c::C,a) = coloralpha(base_color_type(C))(color(c), a)


text:
using KaTeX
using Test

@testset "deps" begin
    @test KaTeX.assets == ["auto-render.min.js", "katex.min.css", "katex.min.js"]
    for asset in KaTeX.assets
        @test isfile(joinpath(KaTeX.assetsdir, asset))
    end
end

const katexfonts = [
    "KaTeX_AMS-Regular",
    "KaTeX_Caligraphic-Bold",
    "KaTeX_Caligraphic-Regular",
    "KaTeX_Fraktur-Bold",
    "KaTeX_Fraktur-Regular",
    "KaTeX_Main-Regular",
    "KaTeX_Main-Bold",
    "KaTeX_Main-Italic",
    "KaTeX_Main-BoldItalic",
    "KaTeX_Math-Italic",
    "KaTeX_Math-BoldItalic",
    "KaTeX_SansSerif-Bold",
    "KaTeX_SansSerif-Italic",
    "KaTeX_SansSerif-Regular",
    "KaTeX_Script-Regular",
    "KaTeX_Size1-Regular",
    "KaTeX_Size2-Regular",
    "KaTeX_Size3-Regular",
    "KaTeX_Size4-Regular",
    "KaTeX_Typewriter-Regular",
]

const exts =["ttf", "woff", "woff2"]

@testset "katexfonts" begin
    @test isdir(KaTeX.fontsdir)
    for font in katexfonts
        for ext in exts
            @test isfile(joinpath(KaTeX.fontsdir, join([font, ext], '.')))
        end
    end
end


text:
import ComradeBase: AbstractModel, IsPrimitive, NotPrimitive, IsAnalytic, NotAnalytic,
                    visanalytic, imanalytic, isprimitive
import ComradeBase: visibility_point,
                    intensitymap, intensitymap!, intensity_point,
                    flux

export visibility, amplitude, closure_phase, logclosure_amplitude, bispectrum,
       visibilities, amplitudes, closure_phases, logclosure_amplitudes, bispectra,
       flux, intensitymap, intensitymap!, PolarizedModel

abstract type AbstractModelImage{M} <: ComradeBase.AbstractModel end

"""
    $(SIGNATURES)
Computes the complex visibility of model `m` at u,v positions `u,v`

If you want to compute the visibilities at a large number of positions
consider using the `visibilities` function which uses MappedArrays to
compute a lazy, no allocation vector.
"""
@inline function visibility(mimg::M, u, v) where {M}
    #first we split based on whether the model is primitive
    _visibility(isprimitive(M), mimg, u, v)
end

@inline function visibility(mimg, uv::ArrayBaselineDatum)
    return visibility(mimg, uv.u, uv.v)
end

"""
    $(SIGNATURES)
Computes the visibility amplitude of model `m` at u,v positions `u,v`

If you want to compute the amplitudes at a large number of positions
consider using the `amplitudes` function which uses MappedArrays to
compute a lazy, no allocation vector.
"""
@inline function amplitude(model, u, v)
    return abs(visibility(model, u, v))
end

"""
    $(SIGNATURES)
Computes the complex bispectrum of model `m` at the uv-triangle
u1,v1 -> u2,v2 -> u3,v3

If you want to compute the bispectrum over a number of triangles
consider using the `bispectra` function which uses MappedArrays to
compute a lazy, no allocation vector.
"""
@inline function bispectrum(model, u1, v1, u2, v2, u3, v3)
    return visibility(model, u1, v1)*visibility(model, u2, v2)*visibility(model, u3, v3)
end

"""
    $(SIGNATURES)
Computes the closure phase of model `m` at the uv-triangle
u1,v1 -> u2,v2 -> u3,v3

If you want to compute closure phases over a number of triangles
consider using the `closure_phases` function which uses MappedArrays to
compute a lazy, no allocation vector.
"""
@inline function closure_phase(model, u1, v1, u2, v2, u3, v3)
    return angle(bispectrum(model, u1, v1, u2, v2, u3, v3))
end

"""
    $(SIGNATURES)
Computes the log-closure amplitude of model `m` at the uv-quadrangle
u1,v1 -> u2,v2 -> u3,v3 -> u4,v3 using the formula

```math
C = \\log\\left|\\frac{V(u1,v1)V(u2,v2)}{V(u3,v3)V(u4,v4)}\\right|
```

If you want to compute log closure amplitudes over a number of triangles
consider using the `logclosure_amplitudes` function which uses MappedArrays to
compute a lazy, no allocation vector.
"""
@inline function logclosure_amplitude(model, u1, v1, u2, v2, u3, v3, u4, v4)
    a1 = amplitude(model, u1, v1)
    a2 = amplitude(model, u2, v2)
    a3 = amplitude(model, u3, v3)
    a4 = amplitude(model, u4, v4)

    return log(a1*a2/(a3*a4))
end

#=
    Welcome to the trait jungle. What is below is
    how we specify how to evaluate the model
=#
@inline function _visibility(::NotPrimitive, m, u, v)
    return visibility_point(m, u, v)
end

@inline function _visibility(::IsPrimitive, m::M, u, v) where {M}
    _visibility_primitive(visanalytic(M), m, u, v)
end

@inline function _visibility_primitive(::IsAnalytic, mimg, u, v)
    return visibility_point(mimg, u, v)
end

@inline function _visibility_primitive(::NotAnalytic, mimg, u, v)
    return mimg.cache.sitp(u, v)
end

@inline function visibilities(m::M, u::AbstractArray, v::AbstractArray) where {M}
    _visibilities(m, u, v)
    #visibilities(visanalytic(M), m, u, v)
end

function visibilities(m, ac::ArrayConfiguration)
    u, v = getuv(ac)
    return visibilities(m, u, v)
end

function _visibilities(m, u::AbstractArray, v::AbstractArray)
    f(x,y) = visibility(m, x, y)
    return mappedarray(f, u, v)
end

"""
    $(SIGNATURES)
Computes the amplitudes of the model `m` at the u,v positions `u`, `v`.

# Notes
If this is a analytic model this is done lazily so the visibilites are only computed
when accessed. Otherwise for numerical model computed with NFFT this is eager.
"""
function amplitudes(m, u::AbstractArray, v::AbstractArray)
    f(x,y) = amplitude(m, x, y)
    return mappedarray(f, u, v)
end

"""
    $(SIGNATURES)
Computes the bispectra of the model `m` at the
triangles (u1,v1), (u2,v2), (u3,v3).

Note this is done lazily so the bispectra is only computed when accessed.
"""
function bispectra(m,
                    u1::AbstractArray, v1::AbstractArray,
                    u2::AbstractArray, v2::AbstractArray,
                    u3::AbstractArray, v3::AbstractArray
                   )
    f(x1,y1,x2,y2,x3,y3) = bispectra(m, x1, y1, x2, y2, x3, y3)
    return mappedarray(f, u1, v1, u2, v2, u3, v3)
end

"""
    $(SIGNATURES)
Computes the closure phases of the model `m` at the
triangles (u1,v1), (u2,v2), (u3,v3).

Note this is done lazily so the closure_phases is only computed when accessed.
"""
function closure_phases(m,
                        u1::AbstractArray, v1::AbstractArray,
                        u2::AbstractArray, v2::AbstractArray,
                        u3::AbstractArray, v3::AbstractArray
                       )
    f(x1,y1,x2,y2,x3,y3) = closure_phase(m, x1, y1, x2, y2, x3, y3)
    return mappedarray(f, u1, v1, u2, v2, u3, v3)
end

"""
    $(SIGNATURES)
Computes the log closure amplitudes of the model `m` at the
quadrangles (u1,v1), (u2,v2), (u3,v3), (u4, v4).

Note this is done lazily so the log closure amplitude is only computed when accessed.
"""
function logclosure_amplitudes(m,
                               u1::AbstractArray, v1::AbstractArray,
                               u2::AbstractArray, v2::AbstractArray,
                               u3::AbstractArray, v3::AbstractArray,
                               u4::AbstractArray, v4::AbstractArray
                              )
    f(x1,y1,x2,y2,x3,y3,x4,y4) = logclosure_amplitude(m, x1, y1, x2, y2, x3, y3, x4, y4)
    return mappedarray(f, u1, v1, u2, v2, u3, v3, u4, v4)
end

function intensitymap!(::NotAnalytic, img::IntensityMap, m, executor=SequentialEx())
    ny, nx = size(img)
    fovx, fovy = fov(img)
    vis = fouriermap(m, fovx, fovy, nx, ny)
    vis = ifftshift(phasedecenter!(vis, fovx, fovy, nx, ny))
    ifft!(vis)
    for I in CartesianIndices(img)
        img[I] = real(vis[I])/(nx*ny)
    end
end

function intensitymap(::NotAnalytic, m, fovx::Real, fovy::Real, nx::Int, ny::Int; pulse=DeltaPulse(), executor=SequentialEx())
    img = IntensityMap(zeros(ny, nx), fovx, fovy, pulse)
    vis = ifftshift(phasedecenter!(fouriermap(m, fovx, fovy, nx, ny), fovx, fovy, nx, ny))
    ifft!(vis)
    for I in CartesianIndices(img)
        img[I] = real(vis[I])/(nx*ny)
    end
    return img
end

include(joinpath(@__DIR__, "modelimage/modelimage.jl"))
include(joinpath(@__DIR__, "modifiers.jl"))
include(joinpath(@__DIR__, "combinators.jl"))
include(joinpath(@__DIR__, "geometric_models.jl"))
include(joinpath(@__DIR__, "polarized.jl"))
include(joinpath(@__DIR__, "test.jl"))