text:
#=
Copyright (c) 2016, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

function pattern_match_call_filter(linfo,f::GlobalRef, id, cond_e, num_cols,table_cols...)
    s = ""
    if f.name!=:__hpat_filter
        return s
    end
    # its an array of array. array[2:end] and table_cols... notation does that
    all_table_cols = table_cols[1]
    out_table_cols = all_table_cols[1:num_cols]
    in_table_cols = all_table_cols[(num_cols + 1):end]

    # For unique counter variables of filter
    unique_id = string(id)
    # assuming that all columns are of same size in a table
    column1_name = ParallelAccelerator.CGen.from_expr(in_table_cols[1],linfo)
    array_length = column1_name*"_array_length_filter" * unique_id
    s *= "int $array_length = " * column1_name * ".ARRAYLEN();\n"
    # Calculate final filtered array length
    write_index = "write_index_filter" * unique_id
    s *= "int $write_index = 1;\n"
    cond_e_arr = ParallelAccelerator.CGen.from_expr(cond_e, linfo)
    s *= "for (int index = 1 ; index < $array_length + 1 ; index++) { \n"
    s *= "if ( $cond_e_arr.ARRAYELEM(index) ){\n"
    # If condition satisfy copy all columns values
    for col_name in in_table_cols
        arr_col_name = ParallelAccelerator.CGen.from_expr(col_name,linfo)
        s *= "$arr_col_name.ARRAYELEM($write_index) =  $arr_col_name.ARRAYELEM(index); \n"
    end
    s *= "$write_index = $write_index + 1;\n"
    s *= "};\n" # if condition
    s *= "};\n" # for loop
    # After filtering we need to change the size of each array
    # And assign to output filter column tables
    for (index, col_name) in enumerate(in_table_cols)
        arr_col_name = ParallelAccelerator.CGen.from_expr(col_name,linfo)
        out_arr_col_name = ParallelAccelerator.CGen.from_expr(out_table_cols[index],linfo)
        s *= "$arr_col_name.dims[0] =  $write_index - 1; \n"
        s *= "$out_arr_col_name = $arr_col_name ; \n"
    end
    return s
end

function pattern_match_call_filter(linfo,f::Any, id, cond_e, num_cols,table_cols...)
    return ""
end

function pattern_match_call_join_seq(linfo, f::GlobalRef, id, table_new_cols_len, table1_cols_len, table2_cols_len, table_columns...)
    s = ""
    if f.name!=:__hpat_join
        return s
    end
    # its an array of array. array[2:end] and table_cols... notation does that
    table_columns = table_columns[1]
    # extract columns of each table
    table_new_cols = table_columns[1:table_new_cols_len]
    table1_cols = table_columns[table_new_cols_len+1:table_new_cols_len+table1_cols_len]
    table2_cols = table_columns[table_new_cols_len+table1_cols_len+1:end]

    # assuming that all columns are of same size in a table
    # Also output table's length would be sum of both table length
    t1c1_length_join = "t1c1_length_join_$id"
    t2c1_length_join = "t2c1_length_join_$id"
    joined_table_length = "joined_table_length_$id"
    t1_c1_join = ParallelAccelerator.CGen.from_expr(table1_cols[1],linfo)
    t2_c1_join = ParallelAccelerator.CGen.from_expr(table2_cols[1],linfo)
    s *= "int $t1c1_length_join = $t1_c1_join.ARRAYLEN() ;\n "
    s *= "int $t2c1_length_join = $t2_c1_join.ARRAYLEN() ;\n "
    s *= "int $joined_table_length = $t2c1_length_join * $t2c1_length_join ;\n "
    # Instantiation of columns for  output table
    for col_name in table_new_cols
        arr_col_name = ParallelAccelerator.CGen.from_expr(col_name,linfo)
        s *= "$arr_col_name = j2c_array<int64_t>::new_j2c_array_1d(NULL, $joined_table_length);\n"
    end
    # Assuming that join is always on the first column of tables
    # Nested for loop implementation of join
    c_cond_sym = "=="
    table_new_counter_join = "table_new_counter_join_$id"
    s *= "int $table_new_counter_join = 1 ; \n"
    s *= "for (int table1_index = 1 ; table1_index < $t1c1_length_join+1 ; table1_index++) { \n"
    s *= "for (int table2_index = 1 ; table2_index < $t2c1_length_join+1 ; table2_index++) { \n"
    s *= "if ( $t1_c1_join.ARRAYELEM(table1_index) $c_cond_sym  $t2_c1_join.ARRAYELEM(table2_index) ){\n"
    count = 0;
    for (index, col_name) in enumerate(table1_cols)
        table_new_col_name = ParallelAccelerator.CGen.from_expr(table_new_cols[index],linfo)
        table1_col_name = ParallelAccelerator.CGen.from_expr(col_name,linfo)
        s *= "$table_new_col_name.ARRAYELEM($table_new_counter_join) = $table1_col_name.ARRAYELEM(table1_index); \n"
        count = count + 1
    end
    for (index, col_name) in enumerate(table2_cols)
        if index == 1
            continue
        end
        table_new_col_name = ParallelAccelerator.CGen.from_expr(table_new_cols[index+count-1],linfo)
        table2_col_name = ParallelAccelerator.CGen.from_expr(col_name,linfo)
        s *= "$table_new_col_name.ARRAYELEM($table_new_counter_join) =  $table2_col_name.ARRAYELEM(table2_index); \n"
    end
    s *= "$table_new_counter_join++;\n"
    s *= "};\n" # join if condition
    s *= "};\n" # inner for loop
    s *= "};\n" # outer for loop
    # Change the size of each output array
    for col_name in table_new_cols
        arr_col_name = ParallelAccelerator.CGen.from_expr(col_name,linfo)
        s *= "$arr_col_name.dims[0] =  $table_new_counter_join - 1; \n"
    end
    # For debugging
    #s *= "for (int i = 1 ; i < $table_new_counter_join ; i++){ std::cout << psale_itemspss_item_sk.ARRAYELEM(i) << std::endl;}\n"
    return s
end

function pattern_match_call_join_seq(linfo, f::Any, id, table_new_len, table1_len, table2_len, table_columns...)
    return ""
end

function pattern_match_call_join(linfo, f::GlobalRef, id, table_new_cols_len, table1_cols_len, table2_cols_len, table_columns...)
    s = ""
    if f.name!=:__hpat_join
        return s
    end
    # TODO remove join random. Use join id/counter in domain pass and pass to this function
    HPAT_path = joinpath(dirname(@__FILE__), "..")
    HPAT_includes = string("\n#include \"", HPAT_path, "/deps/include/hpat_sort.h\"\n")
    ParallelAccelerator.CGen.addCgenUserOptions(ParallelAccelerator.CGen.CgenUserOptions(HPAT_includes,"",""))

    # its an array of array. array[2:end] and table_cols... notation does that
    table_columns = table_columns[1]
    # extract columns of each table
    table_new_cols = table_columns[1:table_new_cols_len]
    table1_cols = table_columns[table_new_cols_len+1:table_new_cols_len+table1_cols_len]
    table2_cols = table_columns[table_new_cols_len+table1_cols_len+1:end]

    s *= "int join_num_pes_$id;\n"
    s *= "MPI_Comm join_comm_$id;\n"

    if haskey(ENV, "ENABLE_GAAS")
        s *= " join_num_pes_$id = " * ENV["ENABLE_GAAS"] * ";\n"
        s *= " join_comm_$id = __hpat_bridge_comm;\n"
    else
        s *= " join_num_pes_$id = __hpat_num_pes ;\n"
        s *= " join_comm_$id = MPI_COMM_WORLD ;\n"
    end

    if haskey(ENV, "ENABLE_GAAS")
        s *= gatherv_part_of_gass(table1_cols, table2_cols, string(id), linfo)
    end

    # Sending counts for both tables
    s *= "int * scount_t1_$id;\n"
    s *= "int * scount_t2_$id;\n"

    s *= "int * scount_t1_tmp_$id;\n"
    s *= "int * scount_t2_tmp_$id;\n"

    s *= "scount_t1_$id = (int*)malloc(sizeof(int)* join_num_pes_$id);\n"
    s *= "memset (scount_t1_$id, 0, sizeof(int)* join_num_pes_$id);\n"

    s *= "scount_t1_tmp_$id = (int*)malloc(sizeof(int)* join_num_pes_$id);\n"
    s *= "memset (scount_t1_tmp_$id, 0, sizeof(int)* join_num_pes_$id);\n"

    s *= "scount_t2_$id = (int*)malloc(sizeof(int)* join_num_pes_$id);\n"
    s *= "memset (scount_t2_$id, 0, sizeof(int)* join_num_pes_$id);\n"

    s *= "scount_t2_tmp_$id = (int*)malloc(sizeof(int)* join_num_pes_$id);\n"
    s *= "memset (scount_t2_tmp_$id, 0, sizeof(int)* join_num_pes_$id);\n"

    # Receiving counts for both tables
    s *= "int rsize_t1_$id = 0;\n"
    s *= "int rsize_t2_$id = 0;\n"

    s *= "int * rcount_t1_$id;\n"
    s *= "int * rcount_t2_$id;\n"
    s *= "rcount_t1_$id = (int*)malloc(sizeof(int)* join_num_pes_$id);\n"
    s *= "rcount_t2_$id = (int*)malloc(sizeof(int)* join_num_pes_$id);\n"

    # Displacement arrays for both tables
    s *= "int * sdis_t1_$id;\n"
    s *= "int * rdis_t1_$id;\n"
    s *= "sdis_t1_$id = (int*)malloc(sizeof(int)* join_num_pes_$id);\n"
    s *= "rdis_t1_$id = (int*)malloc(sizeof(int)* join_num_pes_$id);\n"

    s *= "int *sdis_t2_$id;\n"
    s *= "int *rdis_t2_$id;\n"
    s *= "sdis_t2_$id = (int*)malloc(sizeof(int)* join_num_pes_$id);\n"
    s *= "rdis_t2_$id = (int*)malloc(sizeof(int)* join_num_pes_$id);\n"

    t1c1_length_join = "t1c1_length_join_$id"
    t2c1_length_join = "t2c1_length_join_$id"

    t1_c1_join = ParallelAccelerator.CGen.from_expr(table1_cols[1],linfo)
    t2_c1_join = ParallelAccelerator.CGen.from_expr(table2_cols[1],linfo)
    s *= "int $t1c1_length_join = $t1_c1_join.ARRAYLEN() ;\n "
    s *= "int $t2c1_length_join = $t2_c1_join.ARRAYLEN() ;\n "

    # Starting for table 1
    s *= "for (int i = 1 ; i <  $t1c1_length_join + 1 ; i++){\n"
    s *= "int node_id = $t1_c1_join.ARRAYELEM(i) % join_num_pes_$id ;\n"
    s *= "scount_t1_$id[node_id]++;"
    s *= "}\n"

    s *= "sdis_t1_$id[0]=0;\n"
    s *= "for(int i=1;i < join_num_pes_$id ;i++){\n"
    s *= "sdis_t1_$id[i]=scount_t1_$id[i-1] + sdis_t1_$id[i-1];\n"
    s *= "}\n"

    s *= "MPI_Alltoall(scount_t1_$id,1,MPI_INT,rcount_t1_$id,1,MPI_INT, join_comm_$id);\n"

    # Declaring temporary buffers
    # Assuming that all of them have same length
    # TODO do insertion sort like a combiner in Hadoop
    for (index, col_name) in enumerate(table1_cols)
        table1_col_name = ParallelAccelerator.CGen.from_expr(col_name,linfo)
        j2c_typ = get_j2c_type_from_array(col_name,linfo)
        table1_col_name_tmp = table1_col_name * "_tmp_join_" * string(id)
        s *= "j2c_array< $j2c_typ > $table1_col_name_tmp = j2c_array<$j2c_typ>::new_j2c_array_1d(NULL, $t1c1_length_join );\n"
    end
    s *= "for (int i = 1 ; i <  $t1c1_length_join + 1 ; i++){\n"
    s *= "int node_id = $t1_c1_join.ARRAYELEM(i) % join_num_pes_$id;\n"
    for (index, col_name) in enumerate(table1_cols)
        table1_col_name = ParallelAccelerator.CGen.from_expr(col_name,linfo)
        table1_col_name_tmp = table1_col_name * "_tmp_join_" * string(id)
        s *= "$table1_col_name_tmp.ARRAYELEM(sdis_t1_$id[node_id]+scount_t1_tmp_$id[node_id]+1) = $table1_col_name.ARRAYELEM(i);\n"
    end
    s *= "scount_t1_tmp_$id[node_id]++;\n"
    s *= "}\n"

    # Starting for table 2
    s *= "for (int i = 1 ; i <  $t2c1_length_join + 1 ; i++){\n"
    s *= "int node_id = $t2_c1_join.ARRAYELEM(i) % join_num_pes_$id ;\n"
    s *= "scount_t2_$id[node_id]++;"
    s *= "}\n"

    s *= "sdis_t2_$id[0]=0;\n"
    s *= "for(int i=1;i < __hpat_num_pes;i++){\n"
    s *= "sdis_t2_$id[i]=scount_t2_$id[i-1] + sdis_t2_$id[i-1];\n"
    s *= "}\n"

    s *= "MPI_Alltoall(scount_t2_$id,1,MPI_INT, rcount_t2_$id,1,MPI_INT,MPI_COMM_WORLD);\n"

    # Declaring temporary buffers
    for (index, col_name) in enumerate(table2_cols)
        table2_col_name =ParallelAccelerator.CGen.from_expr(col_name,linfo)
        j2c_typ = get_j2c_type_from_array(col_name,linfo)
        table2_col_name_tmp =  table2_col_name * "_tmp_join_" * string(id)
        s *= "j2c_array< $j2c_typ > $table2_col_name_tmp = j2c_array<$j2c_typ>::new_j2c_array_1d(NULL, $t2c1_length_join);\n"
    end

    s *= "for (int i = 1 ; i <   $t2c1_length_join + 1 ; i++){\n"
    s *= "int node_id = $t2_c1_join.ARRAYELEM(i) % join_num_pes_$id ;\n"
    for (index, col_name) in enumerate(table2_cols)
        table2_col_name =ParallelAccelerator.CGen.from_expr(col_name,linfo)
        table2_col_name_tmp =  table2_col_name * "_tmp_join_$id"
        s *= "$table2_col_name_tmp.ARRAYELEM(sdis_t2_$id[node_id]+scount_t2_tmp_$id[node_id]+1) = $table2_col_name.ARRAYELEM(i);\n"
    end
    s *= "scount_t2_tmp_$id[node_id]++;\n"
    s *= "}\n"

    # Caculating displacements for both tables
    s *= """
              rdis_t1_$id[0]=0;
              rdis_t2_$id[0]=0;
              for(int i=1;i < join_num_pes_$id ;i++){
                  rdis_t1_$id[i] = rcount_t1_$id[i-1] + rdis_t1_$id[i-1];
                  rdis_t2_$id[i] = rcount_t2_$id[i-1] + rdis_t2_$id[i-1];
              }
        """

    # Summing up receiving counts
    s *= """
            for(int i=0;i< join_num_pes_$id ;i++){
                rsize_t1_$id = rsize_t1_$id + rcount_t1_$id[i];
                rsize_t2_$id = rsize_t2_$id + rcount_t2_$id[i];
              }
        """
    for (index, col_name) in enumerate(table1_cols)
        table1_col_name = ParallelAccelerator.CGen.from_expr(col_name,linfo)
        j2c_typ = get_j2c_type_from_array(col_name,linfo)
        table1_col_name_tmp = table1_col_name *"_tmp_join_" * string(id)
        s *= " j2c_array< $j2c_typ > rbuf_$table1_col_name = j2c_array<$j2c_typ>::new_j2c_array_1d(NULL, rsize_t1_$id);\n"
        s *= """ MPI_Alltoallv($table1_col_name_tmp.getData(), scount_t1_$id, sdis_t1_$id, MPI_INT64_T,
                                     rbuf_$table1_col_name.getData(), rcount_t1_$id, rdis_t1_$id, MPI_INT64_T, join_comm_$id);
                     """
        s *= " $table1_col_name = rbuf_$table1_col_name; \n"
    end
    # delete [] tmp_table1_col_name

    for (index, col_name) in enumerate(table2_cols)
        table2_col_name = ParallelAccelerator.CGen.from_expr(col_name,linfo)
        j2c_typ = get_j2c_type_from_array(col_name,linfo)
        table2_col_name_tmp = table2_col_name * "_tmp_join_" * string(id)
        s *= " j2c_array< $j2c_typ > rbuf_$table2_col_name = j2c_array<$j2c_typ>::new_j2c_array_1d(NULL, rsize_t2_$id);\n"
        s *= """ MPI_Alltoallv($table2_col_name_tmp.getData(), scount_t2_$id, sdis_t2_$id, MPI_INT64_T,
                                     rbuf_$table2_col_name.getData(), rcount_t2_$id, rdis_t2_$id, MPI_INT64_T, join_comm_$id);
                     """
        s *= " $table2_col_name = rbuf_$table2_col_name; \n"
    end
    # delete [] tmp_table2_col_name

    table_new_counter_join = "table_new_counter_join_$id"
    s *= "int $table_new_counter_join = 1 ; \n"
    count = 0;
    # Initiatilizing new table(output table) arrays
    for (index, col_name) in enumerate(table1_cols)
        table_new_col_name = ParallelAccelerator.CGen.from_expr(table_new_cols[index],linfo)
        j2c_typ = get_j2c_type_from_array(col_name,linfo)
        s *= "std::vector< $j2c_typ> *vec_$(id)_$table_new_col_name = new std::vector< $j2c_typ>();\n"
        count = count + 1
    end
    for (index, col_name) in enumerate(table2_cols)
        if index == 1
            continue
        end
        table_new_col_name = ParallelAccelerator.CGen.from_expr(table_new_cols[index+count-1],linfo)
        j2c_typ = get_j2c_type_from_array(col_name,linfo)
        s *= "std::vector< $j2c_typ> *vec_$(id)_$table_new_col_name = new std::vector< $j2c_typ>();\n"
    end

    if haskey(ENV, "ENABLE_GAAS")
        s *= scatterv_part_of_gass(table1_cols, table2_cols, string(id), linfo)
    end

    # Use any sorting algorithm here before merging
    # Right now using simple bubble sort
    # TODO add tim sort here too
    j2c_type_t1 = get_j2c_type_from_array(table1_cols[1],linfo)
    j2c_type_t2 = get_j2c_type_from_array(table2_cols[1],linfo)
    t1_length = length(table1_cols)
    t2_length = length(table2_cols)
    t1_all_arrays = "t1_all_arrays_$id"
    t2_all_arrays = "t2_all_arrays_$id"
    s *= "$j2c_type_t1 * $t1_all_arrays[$t1_length - 1];\n"
    s *= "$j2c_type_t2 * $t2_all_arrays[$t2_length - 1];\n"
    for (index, col_name) in enumerate(table1_cols)
        if index == 1
            continue
        end
        arr_index = index - 2
        table1_col_name = ParallelAccelerator.CGen.from_expr(col_name,linfo)
        s *= "$t1_all_arrays [$arr_index] = ( $j2c_type_t1 *) $table1_col_name.getData();\n"
    end

    for (index, col_name) in enumerate(table2_cols)
        if index == 1
            continue
        end
        arr_index = index - 2
        table2_col_name = ParallelAccelerator.CGen.from_expr(col_name,linfo)
        s *= "$t2_all_arrays [$arr_index] = ( $j2c_type_t2 *) $table2_col_name.getData();\n"
    end
    s *= "__hpat_timsort(( $j2c_type_t1 *) $t1_c1_join.getData(), rsize_t1_$id , $t1_all_arrays, $t1_length - 1);\n"
    s *= "__hpat_timsort(( $j2c_type_t2 *) $t2_c1_join.getData(), rsize_t2_$id , $t2_all_arrays, $t2_length - 1);\n"
    # s *= "__hpat_quicksort($t1_all_arrays,$t1_length - 1, ( $j2c_type_t1 *) $t1_c1_join.getData(), 0, $rsize_t1 - 1);\n"
    # s *= "__hpat_quicksort($t2_all_arrays,$t2_length - 1, ( $j2c_type_t2 *) $t2_c1_join.getData(), 0, $rsize_t2 - 1);\n"

    #s *= "qsort($t2_c1_join.getData(),$rsize_t2, sizeof( $j2c_type_t2 ), __hpat_compare_qsort_$j2c_type_t2);\n"
    # after the arrays has been sorted merge them
    # I used algorithm from here www.dcs.ed.ac.uk/home/tz/phd/thesis/node20.htm
    left = "left_join_table_$id"
    right = "right_join_table_$id"
    s *= "int $left = 1;\n"
    s *= "int $right = 1;\n"
    s *= "while ( ($left < rsize_t1_$id + 1) && ($right < rsize_t2_$id + 1) ){\n"
    s *= "if($t1_c1_join.ARRAYELEM($left) == $t2_c1_join.ARRAYELEM($right)){\n"
    count = 0
    for (index, col_name) in enumerate(table1_cols)
        table_new_col_name = ParallelAccelerator.CGen.from_expr(table_new_cols[index],linfo)
        table1_col_name = ParallelAccelerator.CGen.from_expr(col_name,linfo)
        s *= "vec_$(id)_$table_new_col_name->push_back( $table1_col_name.ARRAYELEM($left) );\n"
        count = count + 1
    end
    for (index, col_name) in enumerate(table2_cols)
        if index == 1
            continue
        end
        table_new_col_name = ParallelAccelerator.CGen.from_expr(table_new_cols[index+count-1],linfo)
        table2_col_name = ParallelAccelerator.CGen.from_expr(col_name,linfo)
        s *= "vec_$(id)_$table_new_col_name->push_back( $table2_col_name.ARRAYELEM($right) ); \n"
    end
    s *= "$table_new_counter_join++;\n"

    s *= "int tmp_$left = $left + 1 ;\n"
    s *= "while((tmp_$left < rsize_t1_$id + 1) && ($t1_c1_join.ARRAYELEM(tmp_$left) == $t2_c1_join.ARRAYELEM($right))){\n"
    count = 0
    for (index, col_name) in enumerate(table1_cols)
        table_new_col_name = ParallelAccelerator.CGen.from_expr(table_new_cols[index],linfo)
        table1_col_name = ParallelAccelerator.CGen.from_expr(col_name,linfo)
        s *= "vec_$(id)_$table_new_col_name->push_back( $table1_col_name.ARRAYELEM(tmp_$left) ); \n"
        count = count + 1
    end
    for (index, col_name) in enumerate(table2_cols)
        if index == 1
            continue
        end
        table_new_col_name = ParallelAccelerator.CGen.from_expr(table_new_cols[index+count-1],linfo)
        table2_col_name = ParallelAccelerator.CGen.from_expr(col_name,linfo)
        s *= "vec_$(id)_$table_new_col_name->push_back( $table2_col_name.ARRAYELEM($right) ); \n"
    end
    s *= "tmp_$left++;\n"
    s *= "$table_new_counter_join++;\n"
    s *= "}\n"

    s *= "int tmp_$right = $right + 1 ;\n"
    s *= "while((tmp_$right < rsize_t2_$id + 1) && ($t1_c1_join.ARRAYELEM($left) == $t2_c1_join.ARRAYELEM(tmp_$right))){\n"
    count = 0
    for (index, col_name) in enumerate(table1_cols)
        table_new_col_name = ParallelAccelerator.CGen.from_expr(table_new_cols[index],linfo)
        table1_col_name = ParallelAccelerator.CGen.from_expr(col_name,linfo)
        s *= "vec_$(id)_$table_new_col_name->push_back( $table1_col_name.ARRAYELEM($left) ); \n"
        count = count + 1
    end
    for (index, col_name) in enumerate(table2_cols)
        if index == 1
            continue
        end
        table_new_col_name = ParallelAccelerator.CGen.from_expr(table_new_cols[index+count-1],linfo)
        table2_col_name = ParallelAccelerator.CGen.from_expr(col_name,linfo)
        s *= "vec_$(id)_$table_new_col_name->push_back( $table2_col_name.ARRAYELEM(tmp_$right) ); \n"
    end
    s *= "tmp_$right++;\n"
    s *= "$table_new_counter_join++;\n"
    s *= "}\n"

    s *= "$left++;\n"
    s *= "$right++;\n"
    s *= "}\n" # if condition
    s *= "else if ($t1_c1_join.ARRAYELEM($left) < $t2_c1_join.ARRAYELEM($right))\n"
    s *= "$left++;\n"
    s *= "else\n"
    s *= "$right++;\n"
    s *= "}\n" # while condition

    count = 0
    # Initializing new table(output table) arrays and copy from values from vectors
    for (index, col_name) in enumerate(table1_cols)
        table_new_col_name = ParallelAccelerator.CGen.from_expr(table_new_cols[index],linfo)
        j2c_typ = get_j2c_type_from_array(col_name,linfo)
        s *= "$table_new_col_name = j2c_array<$j2c_typ>::new_j2c_array_1d(vec_$(id)_$table_new_col_name->data(), vec_$(id)_$table_new_col_name->size() );\n"
        count = count + 1
    end
    for (index, col_name) in enumerate(table2_cols)
        if index == 1
            continue
        end
        table_new_col_name = ParallelAccelerator.CGen.from_expr(table_new_cols[index+count-1],linfo)
        j2c_typ = get_j2c_type_from_array(col_name,linfo)
        s *= "$table_new_col_name = j2c_array<$j2c_typ>::new_j2c_array_1d(vec_$(id)_$table_new_col_name->data(), vec_$(id)_$table_new_col_name->size() );\n"
    end
    return s
end

function pattern_match_call_join(linfo, f::Any, id, table_new_len, table1_len, table2_len, table_columns...)
    return ""
end

function pattern_match_call_agg_seq(linfo, f::GlobalRef,  id, groupby_key, num_exprs, expr_func_output_list...)
    s = ""
    if f.name!=:__hpat_aggregate
        return s
    end
    HPAT_path = joinpath(dirname(@__FILE__), "..")
    HPAT_includes = string("\n#include <unordered_map>\n")
    ParallelAccelerator.CGen.addCgenUserOptions(ParallelAccelerator.CGen.CgenUserOptions(HPAT_includes,"",""))

    expr_func_output_list = expr_func_output_list[1]
    exprs_list = expr_func_output_list[1:num_exprs]
    funcs_list = expr_func_output_list[num_exprs+1:(2*num_exprs)]

    # first element of output list is the groupbykey column
    output_cols_list = expr_func_output_list[(2*num_exprs)+1 : end]
    agg_key_col_input = ParallelAccelerator.CGen.from_expr(groupby_key, linfo)
    agg_key_col_output = ParallelAccelerator.CGen.from_expr(output_cols_list[1], linfo)
    # Temporaty map for each column
    for (index, value) in enumerate(output_cols_list)
        table_new_col_name = ParallelAccelerator.CGen.from_expr(value,linfo)
        s *= "std::unordered_map<int,int> temp_map_$table_new_col_name ;\n"
    end
    agg_key_map_temp = "temp_map_$agg_key_col_output"
    s *= "for(int i = 1 ; i < $agg_key_col_input.ARRAYLEN() + 1 ; i++){\n"
    s *= "$agg_key_map_temp[$agg_key_col_input.ARRAYELEM(i)] = $agg_key_col_input.ARRAYELEM(i);\n"
    for (index, func) in enumerate(funcs_list)
        column_name = ""
        expr_name = ParallelAccelerator.CGen.from_expr(exprs_list[index],linfo)
        map_name = "temp_map_" * ParallelAccelerator.CGen.from_expr(output_cols_list[index+1],linfo)
        s *= return_reduction_string_with_closure(agg_key_col_input, expr_name, map_name, func)
    end
    s *= "}\n"
    # Initializing new columns
    for col_name in output_cols_list
        arr_col_name = ParallelAccelerator.CGen.from_expr(col_name, linfo)
        j2c_typ = get_j2c_type_from_array(col_name,linfo)
        s *= "$arr_col_name = j2c_array<$j2c_typ>::new_j2c_array_1d(NULL, $agg_key_map_temp.size());\n"
    end
    # copy back the values from map into arrays
    counter_agg = "counter_agg$id"
    s *= "int $counter_agg = 1;\n"
    s *= "for(auto i : $agg_key_map_temp){\n"
    for (index, value) in enumerate(output_cols_list)
        map_name = ParallelAccelerator.CGen.from_expr(value, linfo)
        s *= "$map_name.ARRAYELEM($counter_agg) = temp_map_$map_name[i.first];\n"
    end
    s *= "$counter_agg++;\n"
    s *= "}\n"
    # Debugging
    # s *= "for (int i = 1 ; i < $counter_agg ; i++){ std::cout << pcustomer_i_classpid3.ARRAYELEM(i) << std::endl;}\n"
    return s
end

function pattern_match_call_agg_seq(linfo, f::Any, id,  groupby_key, num_exprs, exprs_func_list...)
    return ""
end

function pattern_match_call_agg(linfo, f::GlobalRef,  id, groupby_key, num_exprs, expr_func_output_list...)
    s = ""
    if f.name!=:__hpat_aggregate
        return s
    end

    HPAT_path = joinpath(dirname(@__FILE__), "..")
    HPAT_includes = string("\n#include <unordered_map>\n")
    ParallelAccelerator.CGen.addCgenUserOptions(ParallelAccelerator.CGen.CgenUserOptions(HPAT_includes,"",""))

    expr_func_output_list = expr_func_output_list[1]
    exprs_list = expr_func_output_list[1:num_exprs]
    expr_arrs = map(x->ParallelAccelerator.CGen.from_expr(x, linfo), exprs_list)
    funcs_list = expr_func_output_list[num_exprs+1:(2*num_exprs)]

    # first element of output list is the groupbykey column
    output_cols_list = expr_func_output_list[(2*num_exprs)+1 : end]
    agg_key_col_input = ParallelAccelerator.CGen.from_expr(groupby_key, linfo)
    agg_key_col_output = ParallelAccelerator.CGen.from_expr(output_cols_list[1], linfo)

    agg_key_map_temp = "temp_map_$(id)_$agg_key_col_output"
    s *= "std::unordered_map<int,int> $agg_key_map_temp;\n"
    # generate code to get number of elements to send and receive in alltoallv
    # scount_id, rcount_id sdis_id, rdis_id, rsize_id are generated
    s *= gen_alltoall_counts(id, agg_key_col_input, agg_key_map_temp)

    # allocate buffer arrays, first column is groupbykey which is handled separately
    agg_key_col_input_tmp = agg_key_col_input * "_tmp_agg_$id"
    key_ctype = get_j2c_type_from_array(groupby_key, linfo)
    s *= "j2c_array< $key_ctype > $agg_key_col_input_tmp = j2c_array< $key_ctype >::new_j2c_array_1d(NULL, agg_total_unique_keys_$id);\n"
    for (index, expr_arr) in enumerate(expr_arrs)
        expr_arr_tmp = expr_arr * "_tmp_agg_$id"
        j2c_type = get_j2c_type_from_array(output_cols_list[index+1],linfo)
        # allocate aggregate array to accumulate locally
        s *= alloc_agg_arr(funcs_list[index], expr_arr_tmp, key_ctype, j2c_type, "agg_total_unique_keys_$id")
    end

    # aggregate locally and write to send buffer (proper index for each node)
    s *= "for(int i = 1; i <= $agg_key_col_input.ARRAYLEN(); i++){\n"
    s *= "  $key_ctype key = $agg_key_col_input.ARRAYELEM(i);\n"
    s *= "  int node_id = key % __hpat_num_pes ;\n"
    s *= "  if ($agg_key_map_temp.find(key) == $agg_key_map_temp.end()){\n"
    agg_write_index = "agg_write_index_$id"
    s *= "    int $agg_write_index = sdis_$id[node_id]+s_ind_tmp_$id[node_id]+1 ;\n"
    s *= "    $agg_key_map_temp[key] = $agg_write_index;\n"
    s *= "    $agg_key_col_input_tmp.ARRAYELEM($agg_write_index) = key;\n"
    for (index, func) in enumerate(funcs_list)
        expr_arr = expr_arrs[index]
        expr_arr_tmp = expr_arr * "_tmp_agg_$id"
        j2c_type = get_j2c_type_from_array(output_cols_list[index+1],linfo)
        s *= return_combiner_string_with_closure_first_elem(expr_arr_tmp, expr_arr, func, agg_write_index, "key", key_ctype, j2c_type)
    end
    s *= "    s_ind_tmp_$id[node_id]++;\n"
    s *= "}\n"
    s *= "else{\n"
    current_write_index = "current_write_index$id"
    s *= "int $current_write_index = $agg_key_map_temp[key];\n"
    for (index, func) in enumerate(funcs_list)
        expr_arr = expr_arrs[index]
        expr_arr_tmp = expr_arr * "_tmp_agg_$id"
        j2c_type = get_j2c_type_from_array(output_cols_list[index+1],linfo)
        s *= return_combiner_string_with_closure_second_elem(expr_arr_tmp, expr_arr, func, current_write_index, "key", key_ctype, j2c_type)
    end
    s *= "}\n"
    s *= "}\n"
    s *= "$agg_key_map_temp.clear();\n"

    # generate communication for arrays
    # First column is groupbykey which is handled separately
    # After mpi_alltoallv the length of agg_key_col_input is changed. Don't use agg_key_col_input_len
    mpi_type = get_mpi_type_from_array(groupby_key,linfo)
    s *= " j2c_array< $key_ctype > rbuf_$(id)_$agg_key_col_input = j2c_array< $key_ctype >::new_j2c_array_1d(NULL, rsize_$id);\n"
    s *= """ MPI_Alltoallv($agg_key_col_input_tmp.getData(), scount_$id, sdis_$id, $mpi_type,
                                         rbuf_$(id)_$agg_key_col_input.getData(), rcount_$id, rdis_$id, $mpi_type, MPI_COMM_WORLD);
                         """

    for (index, expr_arr) in enumerate(expr_arrs)
        mpi_type = get_mpi_type_from_array(output_cols_list[index+1], linfo)
        j2c_type = get_j2c_type_from_array(output_cols_list[index+1], linfo)
        s *= gen_expr_arr_comm(id, funcs_list[index], expr_arr, j2c_type, mpi_type, key_ctype)
    end
    # delete [] expr_name_tmp

    for col_name in output_cols_list
        j2c_type = get_j2c_type_from_array(col_name,linfo)
        arr_col_name = ParallelAccelerator.CGen.from_expr(col_name, linfo)
        s *= "$arr_col_name = j2c_array< $j2c_type >::new_j2c_array_1d(NULL, rbuf_$(id)_$agg_key_col_input.ARRAYLEN());\n"
    end

    agg_write_index = "agg_write_index_$id"
    s *= "int $agg_write_index = 1;\n"
    s *= "for(int i = 1 ; i < rbuf_$(id)_$agg_key_col_input.ARRAYLEN() + 1 ; i++){\n"
    s *= "  $key_ctype key = rbuf_$(id)_$agg_key_col_input.ARRAYELEM(i);\n"
    s *= "  if ($agg_key_map_temp.find(key) == $agg_key_map_temp.end()){"
    s *= "    $agg_key_map_temp[key] = $agg_write_index ;\n"
    col_name = ParallelAccelerator.CGen.from_expr(output_cols_list[1], linfo)
    s *= "    $col_name.ARRAYELEM($agg_write_index) = key;\n"
    for (index, func) in enumerate(funcs_list)
        expr_arr = expr_arrs[index]
        rbuf_expr_arr = "rbuf_$(id)_" * expr_arr
        new_col_name = ParallelAccelerator.CGen.from_expr(output_cols_list[index + 1], linfo)
        s *= return_reduction_string_with_closure_first_elem(new_col_name, rbuf_expr_arr, func, agg_write_index, "key", expr_arr * "_tmp_agg_$id")
    end
    s *= "    $agg_write_index++;\n"
    s *= "}\n"
    s *= "else{\n"
    current_write_index = "current_write_index$id"
    s *= "int $current_write_index = $agg_key_map_temp[key];\n"
    for (index, func) in enumerate(funcs_list)
        expr_arr = expr_arrs[index]
        rbuf_expr_arr = "rbuf_$(id)_" * expr_arr
        new_col_name = ParallelAccelerator.CGen.from_expr(output_cols_list[index + 1], linfo)
        s *= return_reduction_string_with_closure_second_elem(new_col_name, rbuf_expr_arr, func, current_write_index, "key", expr_arr * "_tmp_agg_$id")
    end
    s *= "}\n"
    s *= "}\n"
    counter_agg = "counter_agg_$id"
    s *= "int $counter_agg = $agg_key_map_temp.size();\n"
    for col_name in output_cols_list
        j2c_type = get_j2c_type_from_array(col_name,linfo)
        arr_col_name = ParallelAccelerator.CGen.from_expr(col_name, linfo)
        s *= "$arr_col_name.dims[0] = $counter_agg;\n"
    end
    return s
end

function pattern_match_call_agg(linfo, f::Any, groupby_key, num_exprs, exprs_func_list...)
    return ""
end

function alloc_agg_arr(func, expr_arr_tmp, key_ctype, arr_ctyp, num_total_keys)
    s = ""
    if func==GlobalRef(HPAT.DomainPass,:length_unique)
        HPAT_includes = """
                            #include <boost/functional/hash.hpp>
                            #include <utility>
                            #include <unordered_set>
                            """
        ParallelAccelerator.CGen.addCgenUserOptions(ParallelAccelerator.CGen.CgenUserOptions(HPAT_includes,"",""))
        # a map for each key
        s = """std::unordered_set<std::pair<$key_ctype,$arr_ctyp>, boost::hash<std::pair<$key_ctype,$arr_ctyp> > > *unique_set_$expr_arr_tmp =
                 new std::unordered_set<std::pair<$key_ctype,$arr_ctyp>, boost::hash<std::pair<$key_ctype,$arr_ctyp> > >[__hpat_num_pes];
               std::unordered_set<std::pair<$key_ctype,$arr_ctyp>, boost::hash<std::pair<$key_ctype,$arr_ctyp> > > recv_unique_set_$expr_arr_tmp;
               std::unordered_map<$key_ctype,$arr_ctyp> unique_map_$expr_arr_tmp;
        """
    else
        s = "j2c_array< $arr_ctyp > $expr_arr_tmp = j2c_array< $arr_ctyp >::new_j2c_array_1d(NULL, $num_total_keys);\n"
    end
    return s
end

# TODO Combine all below five functions into one.
function return_reduction_string_with_closure(agg_key_col_input,expr_arr,agg_map, func::GlobalRef)
    s = ""
    if func.name==:length
        s *= "if ($agg_map.find($agg_key_col_input.ARRAYELEM(i)) == $agg_map.end())\n"
        s *= "$agg_map[$agg_key_col_input.ARRAYELEM(i)] = 1;\n"
        s *= "else \n"
        s *= "$agg_map[$agg_key_col_input.ARRAYELEM(i)] += 1;\n\n"
    elseif func.name==:sum
        s *= "if ($agg_map.find($agg_key_col_input.ARRAYELEM(i)) == $agg_map.end())\n"
        s *= "$agg_map[$agg_key_col_input.ARRAYELEM(i)] = $expr_arr.ARRAYELEM(i) ;\n"
        s *= "else \n"
        s *= "$agg_map[$agg_key_col_input.ARRAYELEM(i)] +=  $expr_arr.ARRAYELEM(i)  ;\n\n"
    elseif func.name==:maximum
        s *= "if ($agg_map.find($agg_key_col_input.ARRAYELEM(i)) == $agg_map.end())){\n"
        s *= "$agg_map[$agg_key_col_input.ARRAYELEM(i)] = $expr_arr.ARRAYELEM(i) ;}\n"
        s *= "else{ \n"
        s *= "if (agg_map_count[$agg_key_col_input.ARRAYELEM(i)] < $expr_arr.ARRAYELEM(i) ) \n"
        s *= "$agg_map[$agg_key_col_input.ARRAYELEM(i)] = $expr_arr.ARRAYELEM(i) ;}\n\n"
    else
        throw("aggregate function not supported in CGen $func")
    end
    return s
end

function return_combiner_string_with_closure_first_elem(new_column_name, expr_arr, func::GlobalRef, write_index, key,key_ctype, j2c_type)
    s = ""
    if func.name==:length
        s *= "$new_column_name.ARRAYELEM($write_index) = 1;\n"
    elseif func.name==:sum
        s *= "$new_column_name.ARRAYELEM($write_index) = $expr_arr.ARRAYELEM(i);\n"
    elseif func.name==:maximum
        s *= "$new_column_name.ARRAYELEM($write_index) = $expr_arr.ARRAYELEM(i);\n"
    elseif func.name==:length_unique
        # s *= "$new_column_name.ARRAYELEM($write_index) = 1;\n"
        s *= "unique_set_$new_column_name[node_id].insert(std::make_pair($key,$expr_arr.ARRAYELEM(i)));\n"
    else
        throw("aggregate function not found $func")
    end
    return s
end

function return_combiner_string_with_closure_second_elem(new_column_name, expr_arr, func::GlobalRef, current_index, key,key_ctype, j2c_type)
    s = ""
    if func.name==:length
        s *= "$new_column_name.ARRAYELEM($current_index) += 1;\n"
    elseif func.name==:sum
        s *= "$new_column_name.ARRAYELEM($current_index) += $expr_arr.ARRAYELEM(i);\n"
    elseif func.name==:maximum
        s *= "if ($new_column_name.ARRAYELEM($current_index) < $expr_arr.ARRAYELEM(i))\n"
        s *= "$new_column_name.ARRAYELEM($current_index) = $expr_arr.ARRAYELEM(i);\n"
    elseif func.name==:length_unique
        #s *= "if(unique_map_$expr_arr[$key].find($expr_arr.ARRAYELEM(i)) == unique_map_$expr_arr[$key].end()){\n"
        #s *= "unique_map_$expr_arr[$key][$expr_arr.ARRAYELEM(i)] = true;\n"
        #s *= "$new_column_name.ARRAYELEM($current_index) += 1;\n"
        #s *= "}\n"
        s *= "unique_set_$new_column_name[node_id].insert(std::make_pair($key,$expr_arr.ARRAYELEM(i)));\n"
    else
        throw("aggregate function not found $func")
    end
    return s
end

# reduction is after alltoallv so data is already partially aggregated
function return_reduction_string_with_closure_first_elem(new_column_name,expr_arr,func::GlobalRef,write_index, key, old_expr_arr)
    s = ""
    if func.name==:length
        s *= "$new_column_name.ARRAYELEM($write_index) = $expr_arr.ARRAYELEM(i);\n"
    elseif func.name==:sum
        s *= "$new_column_name.ARRAYELEM($write_index) = $expr_arr.ARRAYELEM(i);\n"
    elseif func.name==:maximum
        s *= "$new_column_name.ARRAYELEM($write_index) = $expr_arr.ARRAYELEM(i);\n"
    elseif func.name==:length_unique
        #s *= "$new_column_name.ARRAYELEM($write_index) = $expr_arr.ARRAYELEM(i);\n"
        #s *= "unique_map_$old_expr_arr[$key].clear();\n"
        #s *= "unique_map_$old_expr_arr[$key][$expr_arr.ARRAYELEM(i)] = true;\n"
        s *= "$new_column_name.ARRAYELEM($write_index) = unique_map_$old_expr_arr[$key];\n"
    else
        throw("aggregate reduction function not found $func")
    end
    return s
end

function return_reduction_string_with_closure_second_elem(new_column_name,expr_arr,func::GlobalRef, current_index, key, old_expr_arr)
    s = ""
    if func.name==:length
        s *= "$new_column_name.ARRAYELEM($current_index) += $expr_arr.ARRAYELEM(i);\n"
    elseif func.name==:sum
        s *= "$new_column_name.ARRAYELEM($current_index) +=  $expr_arr.ARRAYELEM(i);\n"
    elseif func.name==:maximum
        s *= "if ($new_column_name.ARRAYELEM($current_index) < $expr_arr.ARRAYELEM(i))\n"
        s *= "$new_column_name.ARRAYELEM($current_index) = $expr_arr.ARRAYELEM(i);\n"
    elseif func.name==:length_unique
        #s *= "if(unique_map_$old_expr_arr[$key].find($expr_arr.ARRAYELEM(i)) != unique_map_$old_expr_arr[$key].end()){\n"
        #s *= "unique_map_$old_expr_arr[$key][$expr_arr.ARRAYELEM(i)] = true;\n"
        #s *= "$new_column_name.ARRAYELEM($current_index) += $expr_arr.ARRAYELEM(i);\n"
        #s *= "}\n"
    else
        throw("aggregate reduction function not found $func")
    end
    return s
end

function gen_alltoall_counts(id, agg_key_col_input, agg_key_map_temp)
    s = ""
    # number of elements to send to each node in alltoallv
    s *= "int *scount_$id = (int*)malloc(sizeof(int)*__hpat_num_pes);\n"
    s *= "memset(scount_$id, 0, sizeof(int)*__hpat_num_pes);\n"

    # running write index for send buffer of each node
    s *= "int *s_ind_tmp_$id = (int*)malloc(sizeof(int)*__hpat_num_pes);\n"
    s *= "memset(s_ind_tmp_$id, 0, sizeof(int)*__hpat_num_pes);\n"

    # Receiving counts
    s *= "int rsize_$id = 0;\n"
    s *= "int *rcount_$id = (int*)malloc(sizeof(int)*__hpat_num_pes);\n"

    # send and receive indices for data sent to and received from each node
    s *= "int *sdis_$id = (int*)malloc(sizeof(int)*__hpat_num_pes);\n"
    s *= "int *rdis_$id = (int*)malloc(sizeof(int)*__hpat_num_pes);\n"

    s *= "int agg_total_unique_keys_$id = 0;\n"
    ### Counting displacements for table
    # count unique keys of each node and total unique keys
    s *= "for (int i=1; i <= $agg_key_col_input.ARRAYLEN(); i++){\n"
    s *= "  if ($agg_key_map_temp.find($agg_key_col_input.ARRAYELEM(i)) == $agg_key_map_temp.end()){\n"
    s *= "     $agg_key_map_temp[$agg_key_col_input.ARRAYELEM(i)] = 1;\n"
    s *= "     int node_id = $agg_key_col_input.ARRAYELEM(i) % __hpat_num_pes;\n"
    s *= "     scount_$id[node_id]++;\n"
    s *= "     agg_total_unique_keys_$id++;\n"
    s *= "  }\n"
    s *= "}\n"

    s *= "sdis_$id[0]=0;\n"
    s *= "for(int i=1;i < __hpat_num_pes;i++){\n"
    s *= "sdis_$id[i] = scount_$id[i-1] + sdis_$id[i-1];\n"
    s *= "}\n"

    s *= "MPI_Alltoall(scount_$id,1,MPI_INT,rcount_$id,1,MPI_INT,MPI_COMM_WORLD);\n"

    s *= "$agg_key_map_temp.clear();\n"
    # Caculating displacements
    s *= """
           rdis_$id[0]=0;
           for(int i=1;i < __hpat_num_pes;i++){
               rdis_$id[i] = rcount_$id[i-1] + rdis_$id[i-1];
           }
        """

    # Summing receiving counts
    s *= """
            for(int i=0;i<__hpat_num_pes;i++){
                rsize_$id = rsize_$id + rcount_$id[i];
            }
        """
    return s
end

function gen_expr_arr_comm(id, func, expr_arr, j2c_type, mpi_type, key_ctype)
    s = ""
    if func==GlobalRef(HPAT.DomainPass,:length_unique)
        # sending elements as pairs of (key,value)
        arr_id = "$(id)_$expr_arr"
        sets = "unique_set_$(expr_arr)_tmp_agg_$id"
        # send and recv sizes, in bytes
        s *= "int el_size_$arr_id = sizeof($key_ctype)+sizeof($j2c_type);\n"
        s *= "int total_send_count_$arr_id = 0;\n"
        s *= "int total_recv_count_$arr_id = 0;\n"
        s *= "int *send_sizes_$arr_id = (int*)malloc(sizeof(int)*__hpat_num_pes);\n"
        s *= "memset(send_sizes_$arr_id, 0, sizeof(int)*__hpat_num_pes);\n"
        s *= "int *recv_sizes_$arr_id = (int*)malloc(sizeof(int)*__hpat_num_pes);\n"
        s *= "memset(recv_sizes_$arr_id, 0, sizeof(int)*__hpat_num_pes);\n"
        s *= "int *send_dis_$arr_id = (int*)malloc(sizeof(int)*__hpat_num_pes);\n"
        s *= "int *recv_dis_$arr_id = (int*)malloc(sizeof(int)*__hpat_num_pes);\n"
        s *= """
              for(int i=0; i<__hpat_num_pes; i++) {
                send_sizes_$arr_id[i] = $sets[i].size() * el_size_$arr_id;
                total_send_count_$arr_id += $sets[i].size();
              }
        """
        s *= "MPI_Alltoall(send_sizes_$arr_id, 1, MPI_INT, recv_sizes_$arr_id, 1, MPI_INT, MPI_COMM_WORLD);\n"
        s *= """
              for(int i=0; i<__hpat_num_pes; i++) {
                total_recv_count_$arr_id += recv_sizes_$arr_id[i];
              }
              total_recv_count_$arr_id /= el_size_$arr_id;
        """
        # send and recv buffers
        s *= "char *send_buf_$arr_id = new char[total_send_count_$arr_id * el_size_$arr_id];\n"
        s *= "char *recv_buf_$arr_id = new char[total_recv_count_$arr_id * el_size_$arr_id];\n"
        # pack data in send buffer
        s *= """
             int curr_buff_loc_$arr_id = 0;
             for(int i=0; i<__hpat_num_pes; i++) {
               send_dis_$arr_id[i] = curr_buff_loc_$arr_id;
               for(auto& x: $sets[i]) {
                 $key_ctype *ptr1 = ($key_ctype*) &send_buf_$arr_id[curr_buff_loc_$arr_id];
                 *ptr1 = x.first;
                 curr_buff_loc_$arr_id += sizeof($key_ctype);
                 $j2c_type *ptr2 = ($j2c_type*) &send_buf_$arr_id[curr_buff_loc_$arr_id];
                 *ptr2 = x.second;
                 curr_buff_loc_$arr_id += sizeof($j2c_type);
               }
              }
              recv_dis_$arr_id[0] = 0;
              for(int i=1; i<__hpat_num_pes; i++) {
                recv_dis_$arr_id[i] = recv_sizes_$arr_id[i-1] + recv_dis_$arr_id[i-1];
              }
        """
        s *= """ MPI_Alltoallv(send_buf_$arr_id, send_sizes_$arr_id, send_dis_$arr_id, MPI_CHAR,
                                             recv_buf_$arr_id, recv_sizes_$arr_id, recv_dis_$arr_id, MPI_CHAR, MPI_COMM_WORLD);
                delete[] send_buf_$arr_id;
             """
        s *= """
              curr_buff_loc_$arr_id = 0;
              for(int i=0; i<total_recv_count_$arr_id; i++) {
                  $key_ctype key = *(($key_ctype*)&recv_buf_$arr_id[curr_buff_loc_$arr_id]);
                  curr_buff_loc_$arr_id += sizeof($key_ctype);
                  $j2c_type val = *(($j2c_type*)&recv_buf_$arr_id[curr_buff_loc_$arr_id]);
                  curr_buff_loc_$arr_id += sizeof($j2c_type);
                  recv_unique_set_$(expr_arr)_tmp_agg_$id.insert(std::make_pair(key, val));
              }
              delete[] recv_buf_$arr_id;
              for(auto& x:recv_unique_set_$(expr_arr)_tmp_agg_$id) {
                  if(unique_map_$(expr_arr)_tmp_agg_$id.find(x.first)==unique_map_$(expr_arr)_tmp_agg_$id.end())
                    unique_map_$(expr_arr)_tmp_agg_$id[x.first] = 1;
                  else
                    unique_map_$(expr_arr)_tmp_agg_$id[x.first] += 1;
              }
              recv_unique_set_$(expr_arr)_tmp_agg_$id.clear();
        """
    else
        s *= " j2c_array< $j2c_type > rbuf_$(id)_$expr_arr = j2c_array< $j2c_type >::new_j2c_array_1d(NULL, rsize_$id);\n"
        s *= """ MPI_Alltoallv($(expr_arr)_tmp_agg_$id.getData(), scount_$id, sdis_$id, $mpi_type,
                                             rbuf_$(id)_$expr_arr.getData(), rcount_$id, rdis_$id, $mpi_type, MPI_COMM_WORLD);
             """
    end
end

function pattern_match_call_rebalance(func::GlobalRef, arr::LHSVar, count::LHSVar, linfo)
    s = ""
    if func==GlobalRef(HPAT.API, :__hpat_arr_rebalance)
        typ = ParallelAccelerator.CGen.getSymType(arr, linfo)
        num_dims = ndims(typ)
        typ = eltype(typ)
        c_typ = ParallelAccelerator.CGen.toCtype(typ)
        c_arr = ParallelAccelerator.CGen.from_expr(arr, linfo)
        c_count = ParallelAccelerator.CGen.from_expr(count, linfo)
        mpi_typ = get_mpi_type_from_var_type(typ)
        s *= "int64_t __hpat_old_size_$c_arr = $c_arr.dims[$num_dims-1];\n"
        # get size of each multidim array row (e.g. row size of matrix)
        s *= "int64_t __hpat_row_size_$c_arr = 1;\n"
        for i in 0:num_dims-2
            s *= "__hpat_row_size_$c_arr *= $c_arr.dims[$i];\n"
        end
        # allocate new array
        s *= "$c_typ *__hpat_tmp_$c_arr = new $c_typ[__hpat_row_size_$c_arr*$c_count];\n"
        # copy old data
        s *= "int64_t __hpat_new_data_ind_$c_arr = 0;\n"
        s *= "#define MIN(x, y) (((x) < (y)) ? (x) : (y))\n"
        s *= "for(int64_t i=0; i<MIN(__hpat_old_size_$c_arr, $c_count); i++) {\n"
        s *= "   for(int64_t j=0; j<__hpat_row_size_$c_arr; j++){\n"
        s *= "    __hpat_tmp_$c_arr[__hpat_new_data_ind_$c_arr] = $c_arr.data[__hpat_new_data_ind_$c_arr];\n"
        s *= "    __hpat_new_data_ind_$c_arr++;\n"
        s *= "    }"
        s *= "}\n"
        # my diff, all diffs
        s *= "int64_t _my_diff_$c_arr = __hpat_old_size_$c_arr-$c_count;\n"
        s *= "int64_t *_all_diff_$c_arr = new int64_t[__hpat_num_pes];\n"
        # s *= "printf(\"__hpat_node_id:%d my_size:%d my_count:%d total_size:%d my_diff:%d\\n\", __hpat_node_id, my_size, my_count, total_size, my_diff);"
        s *= "MPI_Allgather(&_my_diff_$c_arr, 1, MPI_LONG_LONG_INT, _all_diff_$c_arr, 1, MPI_LONG_LONG_INT, MPI_COMM_WORLD);\n"
        # printf("__hpat_node_id:%d all_diff[0]:%d all_diff[1]:%d ... all_diff[n-1]:%d\n", __hpat_node_id, all_diff[0], all_diff[1], all_diff[num_pes-1]);
        s *= "MPI_Request *_all_reqs_$c_arr = new MPI_Request[__hpat_num_pes];\n"
        s *= "int _curr_req_$c_arr = 0;\n"
        #// for each potential receiver
        s *= "for(int i=0; i<__hpat_num_pes; i++) {\n"
        #// if receiver
        s *= "  if(_all_diff_$c_arr[i]<0) {\n"
        #// for each potential sender
        s *= "    for(int j=0; j<__hpat_num_pes; j++) {\n"
        #// if sender
        s *= "      if(_all_diff_$c_arr[j]>0) {\n"
        s *= "         int _send_size = MIN(_all_diff_$c_arr[j], -_all_diff_$c_arr[i]);\n"
        #// if I'm receiver
        s *= "         if(__hpat_node_id==i) {\n"
        #//printf("__hpat_node_id:%d receiving from:%d size:%d\n", __hpat_node_id, j, send_size);
        s *= "            MPI_Irecv(&__hpat_tmp_$c_arr[__hpat_new_data_ind_$c_arr], __hpat_row_size_$c_arr*_send_size, $mpi_typ, j, 0, MPI_COMM_WORLD, &_all_reqs_$c_arr[_curr_req_$c_arr++]);\n"
        s *= "            __hpat_new_data_ind_$c_arr += __hpat_row_size_$c_arr*_send_size;\n"
        s *= "         }\n"
        s *= "         if(__hpat_node_id==j) {\n"
        #s *= "            printf("rank:%d sending to:%d size:%d\n", __hpat_node_id, i, send_size);
        s *= "            MPI_Isend(&$c_arr.data[__hpat_new_data_ind_$c_arr], __hpat_row_size_$c_arr*_send_size, $mpi_typ, i, 0, MPI_COMM_WORLD, &_all_reqs_$c_arr[_curr_req_$c_arr++]);\n"
        s *= "            __hpat_new_data_ind_$c_arr += __hpat_row_size_$c_arr*_send_size;\n"
        s *= "         }\n"
        s *= "         _all_diff_$c_arr[i] += _send_size;\n"
        s *= "         _all_diff_$c_arr[j] -= _send_size;\n"
        s *= "         if(_all_diff_$c_arr[i]==0) break;\n"
        s *= "    }\n"
        s *= "   }\n"
        s *= " }\n"
        s *= "}\n"
        s *= "MPI_Waitall(_curr_req_$c_arr, _all_reqs_$c_arr, MPI_STATUSES_IGNORE);\n"
        s *= "delete[] _all_diff_$c_arr;\n"
        s *= "delete[] _all_reqs_$c_arr;\n"
        # delete old array, assign new
        s *= "delete[] $c_arr.data;\n"
        s *= "$c_arr.data = __hpat_tmp_$c_arr;\n"
    end
    return s
end

function pattern_match_call_rebalance(func::ANY, arr::ANY, count::ANY, linfo)
    return ""
end

function gatherv_part_of_gass(table1_cols, table2_cols, id, linfo)
    s = "\n /* Starting Gatherv of gass */ \n"
    table_cols =  [table1_cols; table2_cols]
    for (index, col_name) in enumerate(table_cols)
        table_col_name = ParallelAccelerator.CGen.from_expr(col_name,linfo)
        table_gather_displs =  table_col_name * "_gather_displs_t_" * id
        table_gather_rcounts = table_col_name * "_gather_rcounts_t_" * id
        table_gather_num = table_col_name * "_gather_num_t_" * id
        # initialize all the variable for tables
        s *= "int * $table_gather_displs , * $table_gather_rcounts , $table_gather_num;\n"
        s *= "$table_gather_num = $table_col_name.ARRAYLEN();\n"
        # For node 0
        s *= "$table_gather_rcounts = (int *)malloc(__hpat_num_pes_local*sizeof(int));\n"
        s *= "memset ($table_gather_rcounts, 0, sizeof(int)*__hpat_num_pes_local);\n"
        s *= "MPI_Gather( &$table_gather_num, 1, MPI_INT, $table_gather_rcounts, 1, MPI_INT, 0, __hpat_local_comm);\n"
        # For node 0
        s *= "$table_gather_displs = (int *)malloc(__hpat_num_pes_local*sizeof(int));\n"
        s *= "$table_gather_displs[0] = 0;\n"
        s *= "for (int i=1; i<__hpat_num_pes_local; ++i) {\n"
        s *= "$table_gather_displs[i] = $table_gather_displs[i-1]+ $table_gather_rcounts[i-1];\n"
        s *= "}\n"
        table_gather_rsize = table_col_name * "_gather_rsize_t_" * id
        s *= "int $table_gather_rsize = 0;\n"
        s *= "for(int i=0;i< __hpat_num_pes_local;i++){\n"
    	s *= "$table_gather_rsize = $table_gather_rsize + $table_gather_rcounts[i];\n"
        s *= "}\n"
        # Uncomment for debugging
        # s *= "if (__hpat_node_id_local == 0)\n"
        # s *= "std::cout <<\"table1 $table_col_name: At node 0 After gather recieved = \" << $table_gather_rsize << std::endl;\n"
        table_gather_rbuf = table_col_name * "_gather_rbuf_t_" * id
        s *= "j2c_array< int64_t > $table_gather_rbuf = j2c_array<int64_t>::new_j2c_array_1d(NULL, $table_gather_rsize);\n"
        s *= "MPI_Gatherv($table_col_name.getData(), $table_gather_num, MPI_INT64_T, $table_gather_rbuf.getData(), $table_gather_rcounts, $table_gather_displs, MPI_INT64_T, 0, __hpat_local_comm);\n"
        s *= "$table_col_name = $table_gather_rbuf;\n"
    end
    return s
end

function scatterv_part_of_gass(table1_cols, table2_cols, id, linfo)
    s = "\n /* Starting Scatterv of gass */ \n"
    table_cols =  [table1_cols; table2_cols]
    for (index, col_name) in enumerate(table_cols)
        table_col_name = ParallelAccelerator.CGen.from_expr(col_name,linfo)
        table_scatter_count =  table_col_name * "_scatter_count_t_" * id
        table_scatter_displs = table_col_name * "_scatter_displs_t_" * id
        s *= "int * $table_scatter_count, * $table_scatter_displs ;\n"
        s *= "$table_scatter_count = (int*)malloc(__hpat_num_pes_local * sizeof(int));\n"
        s *= "memset ($table_scatter_count, 0, sizeof(int)*__hpat_num_pes_local);\n"
        s *= "for (int i = 1 ; i <  $table_col_name.ARRAYLEN() + 1 ; i++){\n"
        s *= "int node_id = $table_col_name.ARRAYELEM(i) % __hpat_num_pes_local ;\n"
        s *= "$table_scatter_count[node_id]++;\n"
        s *= "}\n"
        table_scatter_input = table_col_name * "_scatter_input_t_" * id
        s *= "int $table_scatter_input = 0;\n"
        s *= "MPI_Scatter($table_scatter_count, 1, MPI_INT, &$table_scatter_input,1,MPI_INT,0, __hpat_local_comm);\n"
        s *= "$table_scatter_displs = (int *)malloc(__hpat_num_pes_local*sizeof(int));\n"
        s *= "$table_scatter_displs[0] = 0;\n"
        s *= "for (int i=1; i<__hpat_num_pes_local; ++i) {\n"
        s *= "$table_scatter_displs[i] = $table_scatter_displs[i-1]+$table_scatter_count[i-1];\n"
        s *= "}\n"
        table_scatter_rbuf = table_col_name * "_scatter_rbuf_t_" * id
        s *= "j2c_array< int64_t > $table_scatter_rbuf;\n"
        s *= "$table_scatter_rbuf = j2c_array<int64_t>::new_j2c_array_1d(NULL, $table_scatter_input);\n"
        s *= "MPI_Scatterv($table_col_name.getData(), $table_scatter_count , $table_scatter_displs, MPI_INT64_T, $table_scatter_rbuf.getData(), $table_scatter_input, MPI_INT64_T, 0, __hpat_local_comm);\n"
    end
    return s
end


text:
config = MOI.Test.TestConfig()
optimize!(mock) = MOIU.mock_optimize!(mock, [1.0; zeros(5); 6.0; zeros(12); 6.0; zeros(10)],
    (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => zeros(26),
    (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) => zeros(4),
    (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [zeros(7), zeros(7)],
)
for mock in mocks(optimize!)
    Tests.dsos_concave_then_convex_cubic_test(mock, config)
end
optimize!(mock) = MOIU.mock_optimize!(mock, [1.0; zeros(5); 3.0; zeros(3); 3.0; zeros(7); 3.0; zeros(3); 3.0; zeros(5)],
    (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) => zeros(4),
    (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [zeros(7), zeros(7)],
    (MOI.VectorOfVariables, MOI.RotatedSecondOrderCone) => [zeros(3) for i in 1:8],
)
for mock in mocks(optimize!)
    Tests.sdsos_concave_then_convex_cubic_test(mock, config)
end
optimize!(mock) = MOIU.mock_optimize!(mock, [1.0; zeros(5); 6.0; zeros(8); 6.0; zeros(6)],
    (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) => zeros(4),
    (MOI.VectorOfVariables, MOI.PositiveSemidefiniteConeTriangle) => [zeros(6), zeros(6)],
    (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [zeros(7), zeros(7)],
    (MOI.VectorOfVariables, MOI.RotatedSecondOrderCone) => [zeros(3), zeros(3)],
)
for mock in mocks(optimize!)
    Tests.sos_concave_then_convex_cubic_test(mock, config)
end

ext: jl


text:
normpatch(patch) = scaleminmax.(minimum(patch),maximum(patch)).(patch)
viewpatch(patch) =
   [fill(RGB(1,0,0),(1,size(patch,2)+2));
    fill(RGB(1,0,0),(size(patch,1),1)) Gray.(patch) fill(RGB(1,0,0),(size(patch,1),1));
    fill(RGB(1,0,0),(1,size(patch,2)+2))]
viewpatches(patches,zoom=5) = hvcat(size(patches,2),
    permutedims(viewpatch.(repeat.(patches,inner=(zoom,zoom))),[2,1])...)
#;


text:
using Test
using AMR

@testset "ptp" begin
  a = rand(3,3)
  b = rand(3,3)
  c = a + 1.0im * b
  @test ptp(c,1) == ptp(a,1) + 1.0im * ptp(b,1)
end


text:
# Special functions

# See:
#   Martin Maechler (2012) "Accurately Computing log(1 − exp(− |a|))"
#   http://cran.r-project.org/web/packages/Rmpfr/vignettes/log1mexp-note.pdf

# log(1-exp(x))
# NOTE: different than Maechler (2012), no negation inside parantheses
log1mexp(x::Real) = x >= loghalf ? log(-expm1(x)) : log1p(-exp(x))
# log(1+exp(x))
log1pexp(x::Real) = x > zero(x) ? x + log1p(exp(-x)) : log1p(exp(x))
log1pexp(x::Float64) = x <= 18.0 ? log1p(exp(x)) : x <= 33.3 ? x + exp(-x) : x
log1pexp(x::Float32) = x <= 9f0 ? log1p(exp(x)) : x <= 16f0 ? x + exp(-x) : x
log1pexp(x::Integer) = log1pexp(float(x))
# log(exp(x)-1)
logexpm1(x::Real) = x + log1mexp(-x)
logexpm1(x::Float64) = x <= 18.0 ? log(expm1(x)) : x <= 33.3 ? x - exp(-x) : x
logexpm1(x::Float32) = x <= 9f0 ? log(expm1(x)) : x <= 16f0 ? x - exp(-x) : x
logexpm1(x::Integer) = logexpm1(float(x))

# log(1+x^2)
log1psq(x::FloatingPoint) = (ax = abs(x); ax < maxintfloat(x) ? log1p(ax*ax) : 2*log(ax))

φ(z::Real) = exp(-0.5*z*z)/√2π
logφ(z::Real) = -0.5*(z*z + log2π)

Φ(z::Real) = 0.5*erfc(-z/√2)
Φc(z::Real) = 0.5*erfc(z/√2)
logΦ(z::Real) = z < -1.0 ? log(0.5*erfcx(-z/√2)) - 0.5*z*z : log1p(-0.5*erfc(z/√2))
logΦc(z::Real) = z > 1.0 ? log(0.5*erfcx(z/√2)) - 0.5*z*z : log1p(-0.5*erfc(-z/√2))

import Base.Math.@horner

# Rational approximations for the inverse cdf, from:
#   Wichura, M.J. (1988) Algorithm AS 241: The Percentage Points of the Normal Distribution
#   Journal of the Royal Statistical Society. Series C (Applied Statistics), Vol. 37, No. 3, pp. 477-484
for (fn,arg) in ((:Φinv,:p),(:logΦinv,:logp))
    @eval begin
        function $fn($arg::Float32)
            if $(fn == :Φinv)
                q = p - 0.5f0
            else
                q = exp(logp) - 0.5f0
            end
            if abs(q) <= 0.425f0
                r = 0.180625f0 - q*q
                return q * @horner(r,
                                   3.38713_27179f0,
                                   5.04342_71938f1,
                                   1.59291_13202f2,
                                   5.91093_74720f1,
                                   ) /
                @horner(r,
                        1.0f0,
                        1.78951_69469f1,
                        7.87577_57664f1,
                        6.71875_63600f1)
            else
                if $(fn == :Φinv)
                    if p <= 0f0
                        return p == 0f0 ? -inf(Float32) : nan(Float32)
                    elseif p >= 1f0
                        return p == 1f0 ? inf(Float32) : nan(Float32)
                    end
                    r = sqrt(q < 0f0 ? -log(p) : -log1p(-p))
                else
                    if logp == -Inf
                        return -inf(Float32)
                    elseif logp >= 0f0
                        return logp == 0f0 ? inf(Float32) : nan(Float32)
                    end
                    r = sqrt(qf0 < 0 ? -logp : -log1mexp(logp))
                end
                if r < 5.0f0
                    r -= 1.6f0
                    z = @horner(r,
                                1.42343_72777f0,
                                2.75681_53900f0,
                                1.30672_84816f0,
                                1.70238_21103f-1) /
                    @horner(r,
                            1.0f0,
                            7.37001_64250f-1,
                            1.20211_32975f-1)
                else
                    r -= 5.0f0
                    z = @horner(r,
                                6.65790_51150f0,
                                3.08122_63860f0,
                                4.28682_94337f-1,
                                1.73372_03997f-2) /
                    @horner(r,
                            1.0f0,
                            2.41978_94225f-1,
                            1.22582_02635f-2)
                end
                return copysign(z,q)
            end
        end

        function $fn($arg::Real)
            if $(fn == :Φinv)
                q = p - 0.5
            else
                q = exp(logp) - 0.5
            end
            if abs(q) <= 0.425
                r = 0.180625 - q*q
                return q * @horner(r,
                                   3.38713_28727_96366_6080e0,
                                   1.33141_66789_17843_7745e2,
                                   1.97159_09503_06551_4427e3,
                                   1.37316_93765_50946_1125e4,
                                   4.59219_53931_54987_1457e4,
                                   6.72657_70927_00870_0853e4,
                                   3.34305_75583_58812_8105e4,
                                   2.50908_09287_30122_6727e3) /
                @horner(r,
                        1.0,
                        4.23133_30701_60091_1252e1,
                        6.87187_00749_20579_0830e2,
                        5.39419_60214_24751_1077e3,
                        2.12137_94301_58659_5867e4,
                        3.93078_95800_09271_0610e4,
                        2.87290_85735_72194_2674e4,
                        5.22649_52788_52854_5610e3)
            else
                if $(fn == :Φinv)
                    if p <= 0.0
                        return p == 0.0 ? -inf(Float64) : nan(Float64)
                    elseif p >= 1.0
                        return p == 1.0 ? inf(Float64) : nan(Float64)
                    end
                    r = sqrt(q < 0 ? -log(p) : -log1p(-p))
                else
                    if logp == -Inf
                        return -inf(Float64)
                    elseif logp >= 0.0
                        return logp == 0.0 ? inf(Float64) : nan(Float64)
                    end
                    r = sqrt(q < 0 ? -logp : -log1mexp(logp))
                end
                if r < 5.0
                    r -= 1.6
                    z = @horner(r,
                                1.42343_71107_49683_57734e0,
                                4.63033_78461_56545_29590e0,
                                5.76949_72214_60691_40550e0,
                                3.64784_83247_63204_60504e0,
                                1.27045_82524_52368_38258e0,
                                2.41780_72517_74506_11770e-1,
                                2.27238_44989_26918_45833e-2,
                                7.74545_01427_83414_07640e-4) /
                    @horner(r,
                            1.0,
                            2.05319_16266_37758_82187e0,
                            1.67638_48301_83803_84940e0,
                            6.89767_33498_51000_04550e-1,
                            1.48103_97642_74800_74590e-1,
                            1.51986_66563_61645_71966e-2,
                            5.47593_80849_95344_94600e-4,
                            1.05075_00716_44416_84324e-9)
                else
                    r -= 5.0
                    z = @horner(r,
                                6.65790_46435_01103_77720e0,
                                5.46378_49111_64114_36990e0,
                                1.78482_65399_17291_33580e0,
                                2.96560_57182_85048_91230e-1,
                                2.65321_89526_57612_30930e-2,
                                1.24266_09473_88078_43860e-3,
                                2.71155_55687_43487_57815e-5,
                                2.01033_43992_92288_13265e-7) /
                    @horner(r,
                            1.0,
                            5.99832_20655_58879_37690e-1,
                            1.36929_88092_27358_05310e-1,
                            1.48753_61290_85061_48525e-2,
                            7.86869_13114_56132_59100e-4,
                            1.84631_83175_10054_68180e-5,
                            1.42151_17583_16445_88870e-7,
                            2.04426_31033_89939_78564e-15)
                end
                return copysign(z,q)
            end
        end
    end
end

# Multidimensional gamma / partial gamma function
function lpgamma(p::Int, a::Float64)
    res::Float64 = p * (p - 1.0) / 4.0 * log(pi)
    for ii in 1:p
        res += lgamma(a + (1.0 - ii) / 2.0)
    end
    return res
end


text:
push!(LOAD_PATH, "./")
using CAS_QMDP
using ArgParse

# Default:
	# penalty_action=0.02,
	# penalty_closeness=10.0,
	# penalty_nmac=1000.0,
	# penalty_conflict=1.0

function parse_commandline()
	s = ArgParseSettings()
	@add_arg_table s begin
		"--filename"
			help = "Filename of the alpha matrix"
			arg_type = String
			default = "default_filename"
		"--sample"
			help = "If using random sampled noise"
			arg_type = Bool
			default = false
		"--discount"
			help = "Discount factor"
			arg_type = Float64
			default = 0.95
		"--penalty_action"
			help = "Reward function parameter: penalty on action"
			arg_type = Float64
			default = 5.0
		"--penalty_closeness"
			help = "Reward function parameter: penalty on closeness between agents"
			arg_type = Float64
			default = 1.0
		"--penalty_nmac"
			help = "Reward function parameter: penalty on near mid-air collisions"
			arg_type = Float64
			default = 1000.0
		"--penalty_conflict"
			help = "Reward function parameter: penalty on conflict"
			arg_type = Float64
			default = 50.0
	end
	return parse_args(s)
end

function main()
	parsed_args = parse_commandline()

	filename = parsed_args["filename"]
	if_sample = parsed_args["sample"]
	discount = parsed_args["discount"]
	penalty_action = parsed_args["penalty_action"]
	penalty_closeness = parsed_args["penalty_closeness"]
	penalty_nmac = parsed_args["penalty_nmac"]
	penalty_conflict = parsed_args["penalty_conflict"]

	if filename == "default_filename"
		filename = disc * "_maxdist_2km_nmac_150m" *
			"_Discount_" * string(discount) *
			"_ActWeight_" * string(actual_weight)[1 : min(length(string(actual_weight)), 4)] *
			"_IfSample_" * string(if_sample) *
			"_PenAction_" * string(penalty_action) *
			"_PenConflict_" * string(penalty_conflict) *
			"_PenCloseness_" * string(penalty_closeness) *
			"_PenNmac_" * string(penalty_nmac) *
			sigma_str
	end
	# filename = "test_mid2_2"

	open("./VICAS_parameter_sweep_logs.csv", "a") do io
        	@printf(io, "\n")
            @printf(io, "%s,%s,%f,%f,%f,%f,%f,%f,%f",
            	string(Dates.now()), disc, penalty_action, penalty_closeness, penalty_nmac, penalty_conflict,
            	sigma_v, sigma_alert_deg, sigma_coc_deg)
    end

	policy = qmdp(filename,
		discount=discount,
		if_sample=if_sample,
		penalty_action=penalty_action,
		penalty_conflict=penalty_conflict,
		penalty_closeness=penalty_closeness,
		penalty_nmac=penalty_nmac)
	write_alphas(policy, filename)
end

main()

max_stars_repo_name: parachutel/DeepCorrectionCAS


text:
### Main Module for Biochemical Reaction Networks Equilibrate. Exports appropraite functions ###
module BRNEquilibrate

#Fetches required packages
using DiffEqBase
using HomotopyContinuation
using LinearAlgebra
using DynamicPolynomials
using Plots
using ProgressMeter

#Runs the various file contained in the package.
include("equilibrate.jl")
include("bifurcations.jl")
include("plotting.jl")

#Exports the functions that should be avaiable to the user.
export steady_states, solve_poly, bifurcation_diagram, bifurcation_diagram_2d, stability, steady_states_stabiity
export plot, plot!, plot_bifurcation_diagram, plot_bifurcation_diagram!, plot_bifurcation_diagram_2d, plot_bifurcation_diagram_2d!

end

max_stars_repo_name: TorkelE/BRNEquilibrate

max_issues_repo_name: TorkelE/BRNEquilibrate

max_forks_repo_name: TorkelE/BRNEquilibrate


text:
###################################################
## network/tools/graphdep.jl
## find dependence relations between edges
###################################################

"""
	`findNetworkDependence`: express independent edges as sum of dependent edges
	Args:
		n 			: 	the network
		independent : 	list of independent edges
		dependent   : 	list of the remaining (dependent) edges
		numDeps		: 	optional, sparsity parameter
	Returns:
		dependencies:	a m by n matrix, where m is the number of independent edges and n is the total number of edges
						 each column corresponds to the weights of each independent edge in each edge
		edgeMap		:	dictionary, where keys are (orig, dest) pairs and values are the indices in the sorted list of edges
"""
function findNetworkDependence(n::Network, independent::Vector{Edge}, dependent::Vector{Edge}; numDeps::Int = 3)
	# initialize independent edges
	edgeList = collect(edges(n.graph))
	edgeMap = Dict(edgeList[i] => i for i = eachindex(edgeList))
	dep = zeros(length(independent), length(edgeList))
	dependencies = spzeros(length(independent), length(edgeList))
	#special case of all independent edges
	if length(dependent) == 0
		return speye(length(edgeList)), edgeMap
	end
	# build matrix of dependencies
	A = zeros(length(dependent), length(dependent))
	B = zeros(length(dependent), length(independent))
	for (row, edge) in enumerate(dependent)
		nearEdges = findNearEdges(n, edge)
		nNearEdges = length(nearEdges)
		for edge in nearEdges
			if edge in independent
				B[row, findfirst(isequal(edge), independent)] = 1/nNearEdges
			else
				A[row, findfirst(isequal(edge), dependent)] = 1/nNearEdges
			end
		end
	end
	# perform matrix inversion
	Ainv = (eye(size(A)[1]) - A) ^ (-1)
	# create dependency matrix
	depIndices = [edgeMap[edge] for edge in dependent]
	for (row, edge) in enumerate(independent)
		dep[row, edgeMap[edge]] = 1.0
		dep[row, depIndices] = Ainv * B[:, row]
	end
	# sparsify dependency matrix
	for edge in dependent
		dep[:, edgeMap[edge]] = sparsify(dep[:,edgeMap[edge]], numDeps)
	end
	# Multiply and divide by appropriate distances to transform velocity relations to time relations
	for (row, edge) in enumerate(independent)
		dep[row, depIndices] = dep[row, depIndices] .* 1/n.roads[(src(edge), dst(edge))].distance
	end
	for edge in dependent
		dep[:, edgeMap[edge]] = dep[:,edgeMap[edge]] .* n.roads[(src(edge), dst(edge))].distance
	end
	for i in eachindex(independent), j in eachindex(edgeList)
		if dep[i,j] > 0
			dependencies[i,j] = dep[i,j]
		end
	end
	return dependencies, edgeMap
end

"""
	`sparsify` : given a vector of dependences summing to 1, identify main components, reduce all others to 0 and renormalize
	Args:
		dependency 	:	vector summing to 1
		numDeps	  	:	number of components to keep, > 1
		minSize		: 	minimum value that a component must have (after sparsification) to be kept
	Returns:
		newDependency : vector summing to 1, with only numDeps nonzeros
"""
function sparsify(dependency::Vector{Float64}, numDeps::Int, minSize::Float64 = 0.05)
	# find largest components
	p = reverse(sortperm(dependency))
	newDependency = zeros(length(dependency))
	# renormalize
	total = sum(dependency[p[1:numDeps]])
	for i = 1:numDeps
		newDependency[p[i]] = dependency[p[i]]/total
	end
	# if smallest component is too small, re-sparsify with one less component
	if minimum(newDependency[newDependency .> 0]) < minSize
		return sparsify(dependency, numDeps-1, minSize)
	end
	return newDependency
end

"""
	`simplifyPath`: represent path as weighted sum of independent edges
	Args:
		path 		:	path as (Edge, Float) dictionary
		independent : 	list of independent edges
		dependencies:	matrix returned by findNetworkDependence
		edgeMap		:	dictionary, where keys are edges and values are the indices in the list of edges
	Returns:
		newPath		: 	path as (Edge, Float) dictionary, where the keys are the independent edges and the values are their weights in the path
"""
function simplifyPath(path::Dict{Edge, Float64}, independent::Vector{Edge}, dependencies::AbstractArray{Float64,2}, edgeMap::Dict{Edge, Int})
	newPath = Dict{Edge, Float64}()
	dep = zeros(size(dependencies)[1])
	for edge in keys(path)
		dep += path[edge] * dependencies[:,edgeMap[edge]]
	end
	for (i, weight) in enumerate(dep)
		if weight > 0
			newPath[independent[i]] = weight
		end
	end
	return newPath
end

"""
	`evaluateTime` : given time values for independent edges (times), return time of dependent edge (specified by dependency)
	Args:
		n 				: Network object we are working with
		dependencies	: matrix returned by findNetworkDependence
		times 			: array of edge times
		independent 	: list of independent edges
		edgeMap 		: dictionary, where keys are edges and values are the indices in the list of edges
	Returns:
		newTimes 		: array of edge times, where times of independent edges are the same as in input times, and times of dependent edges are calculated from independent edges
"""
function evaluateTimes(n::Network, dependencies::AbstractArray{Float64,2}, times::AbstractArray{Float64,2}, independent::Vector{Edge}, edgeMap::Dict{Edge, Int})
	edgeList = collect(edges(n.graph))
	newTimes = spzeros(Float64, size(times)[1], size(times)[2])
	for (i, edge) in enumerate(edgeList)
		newTimes[src(edge), dst(edge)] = sum([times[src(ind), dst(ind)] * dependencies[idx,i] for (idx, ind) in enumerate(independent)])
	end
	return newTimes
end

"""
	`updateIndependentEdges`	: given set of paths, update independent set of edges
	Args:
		paths		:	vector of vector of paths as edge=> weight dictionaries
		independent :	list of independent edges
		dependent 	:	list of dependent edges
		numEdges	:	number of edges to remove from independent set
		minIndep	:	minimum number of edges that the independent set must always contain
	Returns:
		newIndependent, newDependent : new sets of independent and dependent edges
"""
function updateIndependentEdges(paths::Vector{Vector{Dict{Edge, Float64}}}, independent::Vector{Edge}, dependent::Vector{Edge}, numEdges::Int = 10, minIndep::Int = 100)
	# check if edges can still be removed
	numToRemove = min(numEdges, length(independent) - minIndep)
	if numToRemove == 0
		return independent, dependent
	end
	indices = Dict(independent[i] => i for i=eachindex(independent))
	totalWeight = zeros(length(independent))
	for pathVector in paths, path in pathVector, edge in keys(path)
		totalWeight[indices[edge]] += path[edge]
	end
	# fgind lowest weighted edges
	p = sortperm(totalWeight)
	newIndependent = independent[p[(numToRemove+1):end]]
	newDependent = dependent
	append!(newDependent, independent[p[1:numToRemove]])
	sort!(newDependent)
	return newIndependent, newDependent
end


text:
using ScikitLearn
using CSV # This is a pacakge we use for loading CSV Files.
using DataFrames
using DataFramesMeta
using MLDataPattern
using ProgressMeter
using PyPlot
@sk_import model_selection:train_test_split
@sk_import linear_model:LogisticRegression
@sk_import neural_network:MLPClassifier
@sk_import metrics:accuracy_score
@sk_import metrics:f1_score
@sk_import metrics:classification_report
@sk_import metrics:recall_score
@sk_import metrics:roc_auc_score
@sk_import metrics:confusion_matrix
@sk_import model_selection:cross_val_score
pathtodata = joinpath("julia-scripts", "model-zoo", "covid.csv")
rawdata = CSV.read("covid.csv", DataFrame)

#=================================================================#
#===================== Data Pre-processing =======================#
#=================================================================#

# getting rid of columns we are not interested in
rawdata = rawdata[
    !, Not([:id, :entry_date, :date_symptoms, :date_died, :patient_type, :icu, :sex])
]

# From the data description document we know that 97, 98 or 99 stands for missing data
# so we are getting rid of rows with those values
datawithoutmissing = @linq where(
    :intubed .!= 97,
    :intubed .!= 99,
    :pneumonia .!= 99,
    :pregnancy .!= 97,
    :pregnancy .!= 98,
    :diabetes .!= 98,
    :copd .!= 98,
    :asthma .!= 98,
    :inmsupr .!= 98,
    :hypertension .!= 98,
    :other_disease .!= 98,
    :cardiovascular .!= 98,
    :obesity .!= 98,
    :tobacco .!= 98,
    :renal_chronic .!= 98,
    :contact_other_covid .!= 99,
    :covid_res .!= 3,
)(rawdata)

for col in names(datawithoutmissing)
    replace!(datawithoutmissing[!, col], 2 => 0)
end

pathtocleandata = joinpath("julia-scripts", "model-zoo", "covid_cleaned.csv")
# Saving the cleaned data
# CSV.write(pathtocleandata, datawithoutmissing)

X = convert(Array, datawithoutmissing[!, Not(:covid_res)])
y = convert(Array, datawithoutmissing[!, :covid_res])

X_train, X_test, y_train, y_test = train_test_split(X, y; test_size=0.33, random_state=42)

#=================================================================#
#========================= Model Defenition ======================#
#=================================================================#

# A Simple Logisitic Regression Model with no penalty
simplelogistic = LogisticRegression(; fit_intercept=true, penalty=:none, max_iter=200)
# A Logisitic Regression Model with L2 penalty
logisticL2 = LogisticRegression(; fit_intercept=true, penalty=:l2, max_iter=200)

# A Simple 1 layer Neural Network with 10 nodes
mlp_1layer10nodes = MLPClassifier(; hidden_layer_sizes=10)

# A Simple 1 layer Neural Network with 30 nodes
mlp_1layer30nodes = MLPClassifier(; hidden_layer_sizes=30)

# A 2 layer Neural Network with 30 nodes each
mlp_2layer30nodes = MLPClassifier(; hidden_layer_sizes=(30, 30))

# A 2 layer Neural Network with 100 nodes each
mlp_2layer100nodes = MLPClassifier(; hidden_layer_sizes=(100, 100))

# A 6 layer deep neural network with varying nodes
mlp_6layer = MLPClassifier(; hidden_layer_sizes=(30, 50, 60, 10, 10, 10))

models = [
    simplelogistic,
    logisticL2,
    mlp_1layer10nodes,
    mlp_1layer30nodes,
    mlp_2layer30nodes,
    mlp_2layer100nodes,
    mlp_6layer,
]

function get_prediction(input)
    return predict(simplelogistic, [input])
end

function get_metrics(model, X_train, y_train)
    y_pred = predict(model, X_train)
    acc = accuracy_score(y_train, y_pred)
    f1 = f1_score(y_train, y_pred)
    precision = precision_score(y_train, y_pred)
    recall = recall_score(y_train, y_pred)
    roc = roc_auc_score(y_train, y_pred)
    confusion = confusion_matrix(y_train, y_pred)
    return precision, recall, confusion
end

#=================================================================#
#========================= Model Training ========================#
#=================================================================#
@showprogress for model in models
    fit!(model, X_train, y_train)
end

get_metrics(simplelogistic, X_test, y_test)
y_pred = predict(simplelogistic, X_train)
print(classification_report(y_train, y_pred))

cross_validated = cross_val_score(
    LogisticRegression(; penalty=:none, max_iter=200), X_train, y_train
)

max_stars_repo_name: coinslab/ComputationalCognitiveModeling

max_issues_repo_name: coinslab/ComputationalCognitiveModeling

max_forks_repo_name: coinslab/ComputationalCognitiveModeling


text:
@testset "convert" begin

one = 1
bf = BigFloat(pi)
df = Double64(bf)
sqrt2 = sqrt(2.0)

@test BigFloat(df) == convert(BigFloat, df)
@test Double64(bf) == convert(Double64, bf)

@test Double32(sqrt2) == convert(Double32, sqrt2)
@test Double64(sqrt2) == convert(Double64, sqrt2)

@test one == convert(Int, Double64(1))
@test BigInt(one) == convert(BigInt, Double32(1))
@test Double16(1) == convert(Double16, BigInt(1))

end # convert

@testset "construct $T" for T in (Double16, Double32, Double64)

    i5 = 5
    five = T(5)

    @test Int64(i5) == Int64(five)
    @test Int32(i5) == Int32(five)
    @test Int16(i5) == Int16(five)
    @test Float64(i5) == Float64(five)
    @test Float32(i5) == Float32(five)
    @test Float16(i5) == Float16(five)
    @test five == T(BigFloat(i5))
    @test five == T(BigInt(i5))

end

@testset "maxintfloat $T" for T in (Double16, Double32, Double64)
    @test isinteger(maxintfloat(T))

    # Previous integer is representable, next integer is not
    @test maxintfloat(T) == (maxintfloat(T) - one(T)) + one(T)
    @test maxintfloat(T) != (maxintfloat(T) + one(T)) - one(T)
end

@testset "predicates $T" for T in (Double16, Double32, Double64)

    @test iszero(zero(T))
    @test isone(one(T))
    @test ispositive(one(T))
    @test isnegative(-one(T))
    @test isnonpositive(-one(T))
    @test isnonpositive(zero(T))
    @test isnonnegative(one(T))
    @test isnonnegative(zero(T))

    @test isinf(T(Inf))
    @test isposinf(T(Inf))
    @test isneginf(T(-Inf))
    @test isnan(T(NaN))

    @test isnormal(one(T))
    @test isinteger(one(T))
    @test isfractional(one(T)/2)

end # predicates

@testset "double $T" for T in (Double16, Double32, Double64)

    @test HI(one(T)) == one(T)
    @test LO(one(T)) == zero(T)
    @test HILO(one(T)) == (one(T), zero(T))

    sqrt2  = sqrt(T(2))
    hi, lo = sqrt2.hi, sqrt2.lo
    hilo   = HILO(sqrt2)

    @test HI(sqrt2) == hi
    @test LO(sqrt2) == lo
    @test HILO(sqrt2) == (hi, lo)

    @test HI(hilo) == hi
    @test LO(hilo) == lo
    @test HILO(hilo) == hilo

end # Double

@testset "compare $T" for T in (Double16, Double32, Double64)

    one = 1
    fp1 = T(1)
    two = 2.0
    fp2 = T(2)
    sqrt2  = sqrt(T(2))
    sqrt5  = sqrt(T(5))

    @test isequal(sqrt2, sqrt2)
    @test isless(sqrt2, sqrt5)

    @test sqrt5 != two
    @test two != sqrt5
    @test two == fp2
    @test fp2 == two
    @test two < sqrt5
    @test !(sqrt5 < two)
    @test two <= sqrt5
    @test !(sqrt5 <= two)
    @test !(two > sqrt5)
    @test sqrt5 > two
    @test !(two >= sqrt5)
    @test sqrt5 >= two

    @test sqrt5 != one
    @test one != sqrt5
    @test fp1 == one
    @test one == fp1
    @test one < sqrt5
    @test !(sqrt5 < one)
    @test one <= sqrt5
    @test !(sqrt5 <= one)
    @test !(one > sqrt5)
    @test sqrt5 > one
    @test !(one >= sqrt5)
    @test sqrt5 >= one

    @test   sqrt2 == sqrt2
    @test   sqrt2 >= sqrt2
    @test   sqrt2 <= sqrt2
    @test !(sqrt2 != sqrt2)
    @test !(sqrt2 >  sqrt2)
    @test !(sqrt2 <  sqrt2)

    @test !(sqrt2 == sqrt5)
    @test !(sqrt2 >= sqrt5)
    @test   sqrt2 <= sqrt5
    @test   sqrt2 != sqrt5
    @test !(sqrt2 >  sqrt5)
    @test   sqrt2 <  sqrt5

    @test !(fp2 == sqrt5)
    @test !(fp2 >= sqrt5)
    @test   fp2 <= sqrt5
    @test   fp2 != sqrt5
    @test !(fp2 >  sqrt5)
    @test   fp2 <  sqrt5

    @test !(sqrt5 == fp2)
    @test !(sqrt5 <= fp2)
    @test   sqrt5 >= fp2
    @test   sqrt5 != fp2
    @test !(sqrt5 <  fp2)
    @test   sqrt5 >  fp2

    @test !(two == sqrt5)
    @test !(two >= sqrt5)
    @test   two <= sqrt5
    @test   two != sqrt5
    @test !(two >  sqrt5)
    @test   two <  sqrt5

    @test !(sqrt5 == two)
    @test !(sqrt5 <= two)
    @test   sqrt5 >= two
    @test   sqrt5 != two
    @test !(sqrt5 <  two)
    @test   sqrt5 >  two

end # compare

@testset "eps $T" for T in (Double16, Double32, Double64)
    fp2 = T(2)
    sqrt2  = sqrt(T(2))

    @test eps(fp2) >= DoubleFloats.ulp(fp2)
    @test eps(sqrt2) >= DoubleFloats.ulp(sqrt2)

    @test nextfloat(fp2) > fp2
    @test nextfloat(sqrt2) > sqrt2

    @test prevfloat(fp2) < fp2
    @test prevfloat(sqrt2) < sqrt2

    @test nextfloat(-fp2) > -fp2
    @test nextfloat(-sqrt2) > -sqrt2

    @test prevfloat(-fp2) < -fp2
    @test prevfloat(-sqrt2) < -sqrt2

    @test nextfloat(fp2,2) > fp2
    @test nextfloat(sqrt2,2) > sqrt2

    @test nextfloat(fp2,3) > nextfloat(fp2,2)
    @test prevfloat(fp2,3) < prevfloat(fp2,2)
end


text:
module TestTwoInertiasAndIdealGearWithUnitsAndUncertainties

using TinyModia
using DifferentialEquations
using ModiaPlot
using Measurements
using Unitful

TwoInertiasAndIdealGear = Model(
    J1 = 0.0025u"kg*m^2",
    J2 = (150.0 ± 20.0)u"kg*m^2",
    r  = 105.0,
    tau_max = 1u"N*m",

    phi2 = Var(start = (0.5 ± 0.05)u"rad"),
    w2   = Var(start = 0.0u"rad/s"),
    tau2 = Var(start = 0u"N*m"),

    equations = :[
        tau = if time < 1u"s"; tau_max elseif time < 2u"s"; 0u"N*m" elseif time < 3u"s"; -tau_max else 0u"N*m" end,

        # inertia1
        w1 = der(phi1),
        J1*der(w1) = tau - tau1,

        # gear
        phi1   = r*phi2,
        r*tau1 = tau2,

        # inertia2
        w2 = der(phi2),
        J2*der(w2) = tau2
    ]
)

twoInertiasAndIdealGear = @instantiateModel(TwoInertiasAndIdealGear, FloatType = Measurement{Float64})

simulate!(twoInertiasAndIdealGear, Tsit5(), stopTime = 4.0, logParameters=true, logStates=true)

plot(twoInertiasAndIdealGear, ["phi2", "w2", "der(w2)"])

end


text:
using Dates: unix2datetime, datetime2unix, julian2datetime, datetime2julian
using Mocking: Mocking, @mock

# UTC is an abstract type defined in Dates, for some reason
const utc_tz = FixedTimeZone("UTC")

"""
    DateTime(zdt::ZonedDateTime) -> DateTime

Create a `DateTime` by dropping the associated time zone. Effectively, this new `DateTime`
is implicitly associated with `timezone(zdt)`.

# Example

```jldoctest
julia> zdt = ZonedDateTime(2014, 5, 30, 21, tz"UTC-4")
2014-05-30T21:00:00-04:00

julia> DateTime(zdt)
2014-05-30T21:00:00
```
"""
Dates.DateTime(zdt::ZonedDateTime) = zdt.utc_datetime + zdt.zone.offset

"""
    DateTime(zdt::ZonedDateTime, ::Type{UTC}) -> DateTime

Create a `DateTime` which is implicitly associated with UTC.

# Example

```jldoctest
julia> zdt = ZonedDateTime(2014, 5, 30, 21, tz"UTC-4")
2014-05-30T21:00:00-04:00

julia> DateTime(zdt, UTC)
2014-05-31T01:00:00
```
"""
Dates.DateTime(zdt::ZonedDateTime, ::Type{UTC}) = zdt.utc_datetime

"""
    Date(zdt::ZonedDateTime) -> Date

Create a `Date` by dropping the associated time zone. Effectively, this new `Date`
is implicitly associated with `timezone(zdt)`.

# Example

```jldoctest
julia> zdt = ZonedDateTime(2014, 5, 30, 21, tz"UTC-4")
2014-05-30T21:00:00-04:00

julia> Date(zdt)
2014-05-30
```
"""
Dates.Date(zdt::ZonedDateTime) = Date(DateTime(zdt))

"""
    Date(zdt::ZonedDateTime, ::Type{UTC}) -> Date

Create a `Date` which is implicitly associated with UTC.

# Example

```jldoctest
julia> zdt = ZonedDateTime(2014, 5, 30, 21, tz"UTC-4")
2014-05-30T21:00:00-04:00

julia> Date(zdt, UTC)
2014-05-31
```
"""
Dates.Date(zdt::ZonedDateTime, ::Type{UTC}) = Date(DateTime(zdt, UTC))

"""
    Time(zdt::ZonedDateTime) -> Time

Create a `Time` by dropping the associated time zone. Effectively, this new `Time`
is implicitly associated with `timezone(zdt)`.

# Example

```jldoctest
julia> zdt = ZonedDateTime(2014, 5, 30, 21, tz"UTC-4")
2014-05-30T21:00:00-04:00

julia> Time(zdt)
21:00:00
```
"""
Dates.Time(zdt::ZonedDateTime) = Time(DateTime(zdt))

"""
    Time(zdt::ZonedDateTime, ::Type{UTC}) -> Date

Create a `Time` which is implicitly associated with UTC.

# Example

```jldoctest
julia> zdt = ZonedDateTime(2014, 5, 30, 21, tz"UTC-4")
2014-05-30T21:00:00-04:00

julia> Time(zdt, UTC)
01:00:00
```
"""
Dates.Time(zdt::ZonedDateTime, ::Type{UTC}) = Time(DateTime(zdt, UTC))

"""
    now(::TimeZone) -> ZonedDateTime

Returns a `ZonedDateTime` corresponding to the user's system time in the specified `TimeZone`.
"""
function Dates.now(tz::TimeZone)
    utc = unix2datetime(time())
    ZonedDateTime(utc, tz, from_utc=true)
end

"""
    today(tz::TimeZone) -> Date

Returns the date portion of `now(tz)` in local time.

# Examples

```julia
julia> a, b = now(tz"Pacific/Midway"), now(tz"Pacific/Apia")
(2017-11-09T03:47:04.226-11:00, 2017-11-10T04:47:04.226+14:00)

julia> a - b
0 milliseconds

julia> today(tz"Pacific/Midway"), today(tz"Pacific/Apia")
(2017-11-09, 2017-11-10)
```
"""
Dates.today(tz::TimeZone) = Date(now(tz))

"""
    todayat(tod::Time, tz::TimeZone, [amb]) -> ZonedDateTime

Creates a `ZonedDateTime` for today at the specified time of day. If the result is ambiguous
in the given `TimeZone` then `amb` can be supplied to resolve ambiguity.

# Examples

```julia
julia> today(tz"Europe/Warsaw")
2017-11-09

julia> todayat(Time(10, 30), tz"Europe/Warsaw")
2017-11-09T10:30:00+01:00
```
"""
function todayat(tod::Time, tz::VariableTimeZone, amb::Union{Integer,Bool})
    ZonedDateTime((@mock today(tz)) + tod, tz, amb)
end

todayat(tod::Time, tz::TimeZone) = ZonedDateTime((@mock today(tz)) + tod, tz)

"""
    astimezone(zdt::ZonedDateTime, tz::TimeZone) -> ZonedDateTime

Converts a `ZonedDateTime` from its current `TimeZone` into the specified `TimeZone`.
"""
function astimezone end

function astimezone(zdt::ZonedDateTime, tz::VariableTimeZone)
    i = searchsortedlast(
        tz.transitions, zdt.utc_datetime,
        by=v -> typeof(v) == Transition ? v.utc_datetime : v,
    )

    if i == 0
        throw(NonExistentTimeError(DateTime(zdt), tz))
    end

    zone = tz.transitions[i].zone
    return ZonedDateTime(zdt.utc_datetime, tz, zone)
end

function astimezone(zdt::ZonedDateTime, tz::FixedTimeZone)
    return ZonedDateTime(zdt.utc_datetime, tz, tz)
end

function zdt2julian(zdt::ZonedDateTime)
    datetime2julian(DateTime(zdt, UTC))
end

function zdt2julian(::Type{T}, zdt::ZonedDateTime) where T<:Integer
    floor(T, datetime2julian(DateTime(zdt, UTC)))
end

function zdt2julian(::Type{T}, zdt::ZonedDateTime) where T<:Real
    convert(T, datetime2julian(DateTime(zdt, UTC)))
end

function julian2zdt(jd::Real)
    ZonedDateTime(julian2datetime(jd), utc_tz, from_utc=true)
end

function zdt2unix(zdt::ZonedDateTime)
    datetime2unix(DateTime(zdt, UTC))
end

function zdt2unix(::Type{T}, zdt::ZonedDateTime) where T<:Integer
    floor(T, datetime2unix(DateTime(zdt, UTC)))
end

function zdt2unix(::Type{T}, zdt::ZonedDateTime) where T<:Real
    convert(T, datetime2unix(DateTime(zdt, UTC)))
end

function unix2zdt(seconds::Real)
    ZonedDateTime(unix2datetime(seconds), utc_tz, from_utc=true)
end


text:
module GraphletCounting

using LightGraphs

include("DegreeOrientedDAG.jl")
export sort_by_deg_or!, toDAG

include("CountTriangles.jl")
export enumerate_wedges

end