/*****************************************************************************
*
* Copyright (c) 2000 - 2017, Lawrence Livermore National Security, LLC
* Produced at the Lawrence Livermore National Laboratory
* LLNL-CODE-442911
* All rights reserved.
*
* This file is  part of VisIt. For  details, see https://visit.llnl.gov/.  The
* full copyright notice is contained in the file COPYRIGHT located at the root
* of the VisIt distribution or at http://www.llnl.gov/visit/copyright.html.
*
* 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 disclaimer below.
*  - Redistributions in binary form must reproduce the above copyright notice,
*    this  list of  conditions  and  the  disclaimer (as noted below)  in  the
*    documentation and/or other materials provided with the distribution.
*  - Neither the name of  the LLNS/LLNL nor the names of  its contributors may
*    be used to endorse or promote products derived from this software without
*    specific prior written permission.
*
* 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 LAWRENCE  LIVERMORE NATIONAL  SECURITY,
* LLC, THE  U.S.  DEPARTMENT OF  ENERGY  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.
*
*****************************************************************************/

// ************************************************************************* //
//                        avtCurveConstructorFilter.C                        //
// ************************************************************************* //

#include <avtCurveConstructorFilter.h>
#include <avtDataTree.h>

#include <vtkCellArray.h>
#include <vtkDataArray.h>
#include <vtkDataSet.h>
#include <vtkPointData.h>
#include <vtkPolyData.h>
#include <vtkRectilinearGrid.h>
#include <vtkVisItUtility.h>
#include <InvalidDimensionsException.h>
#include <NoInputException.h>
#include <MapNode.h>

#include <DebugStream.h>

#ifdef PARALLEL
#include <mpi.h>
#include <avtParallel.h>
#include <vtkDataSetReader.h>
#include <vtkDataSetWriter.h>
#include <vtkCharArray.h>
#endif


#include <maptypes.h>



// ****************************************************************************
//  Method: avtCurveConstructorFilter constructor
//
//  Programmer: kbonnell -- generated by xml2info
//  Creation:   Sat Apr 20 13:01:58 PST 2002
//
//  Modifications:
//    Kathleen Bonnell, Thu Mar 19 17:42:48 PDT 2009
//    Added forceConstruct.
//
//    Kathleen Bonnell, Mon Mar 23 09:53:41 PDT 2009
//    Removed forceConstruct.
//
// ****************************************************************************

avtCurveConstructorFilter::avtCurveConstructorFilter()
{
}


// ****************************************************************************
//  Method: avtCurveConstructorFilter destructor
//
//  Programmer: kbonnell -- generated by xml2info
//  Creation:   Sat Apr 20 13:01:58 PST 2002
//
// ****************************************************************************

avtCurveConstructorFilter::~avtCurveConstructorFilter()
{
}


// ****************************************************************************
//  Method: avtCurveConstructorFilter::ExecuteData
//
//  Purpose:
//      Does the actual VTK code to modify the dataset.
//
//  Arguments:
//      inDS      The input dataset.
//      <unused>  The domain number.
//
//  Returns:      The output dataset.
//
//  Programmer: kbonnell -- generated by xml2info
//  Creation:   Sat Apr 20 13:01:58 PST 2002
//
//  Modifications:
//
//    Hank Childs, Tue May 28 11:41:21 PDT 2002
//    Use the variable name as the label for our output data tree.
//
//    Kathleen Bonnell, Fri Jul 12 16:53:11 PDT 2002
//    Removed vtk filters associated with label-creation.  Now handled by
//    the plot.
//
//    Kathleen Bonnell, Tue Dec 23 10:18:06 PST 2003 
//    Added logic to handle point-sorting when necessary. Add vertex cells,
//    so they can be displayed upon user request.
//    
//    Mark C. Miller, Wed Jun  9 21:50:12 PDT 2004
//    Eliminated use of MPI_ANY_TAG and modified to use GetUniqueMessageTags
//
//    Hank Childs, Wed Jan  4 11:21:59 PST 2006
//    Allocate the size of the verts array front-end.
//
//    Kathleen Bonnell, Tue Jun 20 16:02:38 PDT 2006
//    Save the curve points in output array to be added to PlotInfoAttributes.
//
//    Kathleen Bonnell, Mon Jul 31 16:50:55 PDT 2006 
//    Changed reader from PolyData to DataSetReader, as curves can now be
//    represented as 1D RectilinearGrids.  Modified logic to handle
//    Rectilinear Grids as needed.
//
//    Hank Childs, Thu Jan  4 09:29:59 PST 2007
//    Add handling for error case where a tree is not empty, but it has no
//    leaves (so it really is "IsEmpty").
//
//    Mark C. Miller, Mon Jan 22 22:09:01 PST 2007
//    Changed MPI_COMM_WORLD to VISIT_MPI_COMM
//
//    Dave Bremer, Thu Jun  7 19:47:35 PDT 2007
//    Added a fix to correctly construct a curve from a histogram plot.
//    I determine if there's a problem by looking at the geometry, and
//    if necessary flip the last two points in each group of 4. 
//
//    Hank Childs, Fri Feb 15 15:52:46 PST 2008
//    Throw an exception in an error condition.
//
//    Mark C. Miller, Wed Jun 11 12:09:30 PDT 2008
//    Replaced Exception with continue
//
//    Kathleen Bonnell, Thu Mar 19 17:42:48 PDT 2009
//    Don't use unique-leaf test for early termination if forceConstruct is set.
//
//    Kathleen Bonnell, Mon Mar 23 09:54:09 PDT 2009
//    Removed unique-leaf test and forceConstruct.
//
//    Brad Whitlock, Tue Mar 30 15:33:19 PDT 2010
//    Pass along avtCurveTransform in the field data if it exists.
//
//    Kathleen Bonnell, Tue Dec 14 12:57:14 PST 2010
//    Ensure the output variable is named.
//
//    Kathleen Bonnell, Thu Feb 17 09:19:06 PST 2011
//    Moved bulk of dataset construction into CreateSingleOutput method,
//    which will allow creation of multiple outputs when necessary.
//
//    Kathleen Biagas, Thu Sep 29 06:11:38 PDT 2011
//    Only construct multiple outputs when requested in data attributes.
//
// ****************************************************************************

void avtCurveConstructorFilter::Execute()
{
    avtDataTree_p inTree = GetInputDataTree();

#ifdef PARALLEL
    //
    //  Gather all data onto one processor, so that 
    //  there will be no discontinuities in the curve.
    //
    int myRank, numProcs;

    MPI_Comm_size(VISIT_MPI_COMM, &numProcs);
    MPI_Comm_rank(VISIT_MPI_COMM, &myRank);

    int tags[3];
    GetUniqueMessageTags(tags, 3);
    int mpiNdsTag  = tags[0];
    int mpiSizeTag = tags[1];
    int mpiDataTag = tags[2];

    if (myRank == 0)
    {
        int i , j;
        for (i = 1; i < numProcs; i++)
        {
           MPI_Status stat;
           MPI_Status stat2;
           int nds = 0, size = 0;
           MPI_Recv(&nds, 1, MPI_INT, MPI_ANY_SOURCE, mpiNdsTag,
                    VISIT_MPI_COMM, &stat);
           for (j = 0; j < nds; j++)
           {
               vtkDataSetReader *reader = vtkDataSetReader::New(); 
               reader->ReadFromInputStringOn();
               MPI_Recv(&size, 1, MPI_INT, stat.MPI_SOURCE, mpiSizeTag,
                         VISIT_MPI_COMM, &stat2);
               char *str = new char[size];
               MPI_Recv(str, size, MPI_CHAR, stat.MPI_SOURCE, mpiDataTag,
                        VISIT_MPI_COMM, &stat2);
               vtkCharArray *charArray = vtkCharArray::New();
               charArray->SetArray((char*)str, size, 1);
               reader->SetInputArray(charArray);
               reader->Update();
      
               inTree->Merge(new avtDataTree(reader->GetOutput(), -1));
               delete [] str;
               reader->Delete(); 
               charArray->Delete(); 
           }
        }
    }
    else
    {
        SetOutputDataTree(new avtDataTree());
        char *str = NULL;;
        int i = 0, size = 0, nleaves = 0;
        if (inTree->IsEmpty())
        {
            MPI_Send(&nleaves, 1, MPI_INT, 0, mpiNdsTag, VISIT_MPI_COMM);
            return;
        }

        vtkDataSet **ds = inTree->GetAllLeaves(nleaves);
        MPI_Send(&nleaves, 1, MPI_INT, 0, mpiNdsTag, VISIT_MPI_COMM);
 
        vtkDataSetWriter *writer = vtkDataSetWriter::New();
        writer->WriteToOutputStringOn();
        writer->SetFileTypeToBinary();

        for (i = 0; i < nleaves; i++)
        {
            writer->SetInputData(ds[i]);
            writer->Write();
            size =  writer->GetOutputStringLength();
            str  =  writer->RegisterAndGetOutputString();

            MPI_Send(&size, 1, MPI_INT, 0, mpiSizeTag, VISIT_MPI_COMM);
            MPI_Send(str, size, MPI_CHAR, 0, mpiDataTag, VISIT_MPI_COMM);
            delete [] str; //allocated by writer
        }
        writer->Delete(); 
        delete [] ds;
        return;
    }
#endif

    if (inTree->IsEmpty())
    {
        SetOutputDataTree(inTree);
        return; 
    }

    avtDataTree_p outTree;
    stringVector labels;
    inTree->GetAllLabels(labels);
    
    if (labels.size() == 0 || 
        !GetInput()->GetInfo().GetAttributes().GetConstructMultipleCurves())
    {
        vtkDataSet *outGrid;
        outGrid = CreateSingleOutput(inTree);
        if (outGrid == NULL)
        {
            SetOutputDataTree(inTree);
            return; 
        }
        const char *vname = (pipelineVariable != NULL ? pipelineVariable : "");
        avtDataRepresentation dr(outGrid, -1, vname);
        outTree = new avtDataTree(dr);
        outGrid->Delete();
    }
    else
    {
        vtkDataSet **ds = new vtkDataSet *[labels.size()];
        for (size_t i = 0; i < labels.size(); ++i)
        {
            ds[i] = NULL;
            avtDataTree_p oneTree = inTree->PruneTree(labels[i]);
            ds[i] = CreateSingleOutput(oneTree);
            if (ds[i] == NULL)
            {
                for (size_t j = 0; j < i; ++j)
                {
                    if (ds[j] != NULL)
                        ds[j]->Delete(); 
                }
                delete [] ds;
                SetOutputDataTree(inTree);
                return; 
            }
        }
        outTree = new avtDataTree((int)labels.size(), ds, -1, labels);
    }

    SetOutputDataTree(outTree);
}
 
 
// ****************************************************************************
//  Method: avtCurveConstructorFilter::CreateSingleOutput
//
//  Notes:  Moved from Execute mthod.  
//
//  Purpose:
//      Does the actual VTK code to modify the dataset.
//
//  Arguments:
//      inTree    The input data tree.
//
//  Returns:      The output dataset.
//
//  Programmer: Kathleen Bonnell
//  Creation:   February 17, 2011
//
//  Modifications:
//
// ****************************************************************************

vtkDataSet *
avtCurveConstructorFilter::CreateSingleOutput(avtDataTree_p inTree)
{
    //
    //  This filter doesn't do much right now.  Basically just a 
    //  "connect-the-dots" between the vertices.  
    //
    int nleaves, j, k;
    double x;
    vtkDataSet **ds;
    ds = inTree->GetAllLeaves(nleaves);

    if (nleaves == 0)
    {
        // Kind of a bizarre error, since "IsEmpty" above should have returned
        // true.  This situation does occur, and it is a little dangerous to
        // change the behavior of IsEmpty, so just accomodate.
        delete [] ds;
        return NULL;
    }
    DoubleIntMap minX;

    //
    //  Cannot assume that there are no overlaps of points
    //  between datasets, so keep track of first and last
    //  x values, to determine if sorting is required. 
    //  Now must order the datasets so that the points will
    //  be accessed in proper order, to avoid discontinuities
    //  between domains.
    //  Store in map with x value as the key.  
    //  (from first point in each ds).
    //
    double *mm = new double[nleaves*2];
    int npts;

    for (j = 0; j < nleaves; j++)
    {
        vtkDataArray *data = NULL;
        if (ds[j]->GetDataObjectType() == VTK_RECTILINEAR_GRID)
        {
            data = ((vtkRectilinearGrid*)ds[j])->GetXCoordinates();
        }
        else if (ds[j]->GetDataObjectType() == VTK_POLY_DATA)
        {
            data = ((vtkPolyData*)ds[j])->GetPoints()->GetData();
        }
        else
        {
            continue;
        }
        npts = data->GetNumberOfTuples();
        x = data->GetComponent(0, 0);
        minX.insert(DoubleIntMap::value_type(x, j));
        mm[j*2] = x;
        mm[j*2+1] = data->GetComponent(npts-1, 0); 
    }

    bool requiresSort = false;
    for (j = 0; j < nleaves && !requiresSort; j++)
    {
        for (k = 0; k < nleaves && !requiresSort; k++)
        {
            if (k != j && mm[k*2] > mm[j*2] && mm[k*2] < mm[j*2+1])
                requiresSort = true;
        }
    }
    delete [] mm;

    vtkDataArray *inXC; 
    vtkDataArray *inVal; 

    int dtype = VTK_FLOAT;
    if (ds[0]->GetDataObjectType() == VTK_RECTILINEAR_GRID)
        dtype = ((vtkRectilinearGrid*)ds[0])->GetXCoordinates()->GetDataType();
    else 
        dtype = ((vtkPolyData*)ds[0])->GetPoints()->GetData()->GetDataType();
    vtkRectilinearGrid *outGrid = vtkVisItUtility::Create1DRGrid(0, dtype);
    vtkDataArray *outXC  = outGrid->GetXCoordinates();
    vtkDataArray *outVal = outXC->NewInstance();
    outGrid->GetPointData()->SetScalars(outVal);
    outVal->Delete();
    int nPoints; 

    //
    //  Ensure that the output is 2d by setting z-component to zero.
    //
    DoubleIntMap::iterator it;
    int index = 0;
    for (it = minX.begin(); it != minX.end(); it++)
    {
        if (ds[(*it).second]->GetDataObjectType() == VTK_RECTILINEAR_GRID)
        {
            inXC  = ((vtkRectilinearGrid*)ds[(*it).second])->GetXCoordinates();
            inVal = ((vtkRectilinearGrid*)ds[(*it).second])->GetPointData()->
                    GetScalars();
            nPoints = inXC->GetNumberOfTuples();
            for (vtkIdType i = 0; i < nPoints; i++, index++)
            {
                outXC->InsertNextTuple1(inXC->GetTuple1(i));
                outVal->InsertNextTuple1(inVal->GetTuple1(i));
            }
        }
        else if (ds[(*it).second]->GetDataObjectType() == VTK_POLY_DATA)
        {
            inXC  = ((vtkPolyData*)ds[(*it).second])->GetPoints()->GetData();
            nPoints = inXC->GetNumberOfTuples();
            
            // Insert a check to see if this data comes from a histogram plot,
            // in which case we need to swap some of the points.  Data comes
            // in groups of 4 points for each bar, connected like this: N
            // so we reverse the last two points to create bar shapes instead.
            bool bIsHistogramPlot = false;
            if (nPoints%4 == 0)
            {
                bIsHistogramPlot = true;
                for (int i = 0; i < nPoints; i+=4)
                {
                    if ( inXC->GetComponent(i,   0) != inXC->GetComponent(i+1, 0) ||
                         inXC->GetComponent(i+2, 0) != inXC->GetComponent(i+3, 0) ||
                         inXC->GetComponent(i,   1) != inXC->GetComponent(i+2, 1) ||
                         inXC->GetComponent(i+1, 1) != inXC->GetComponent(i+3, 1) )
                    {
                        bIsHistogramPlot = false;
                        break;
                    }
                }
            }
            if (!bIsHistogramPlot)
            {
                for (int i = 0; i < nPoints; i++, index++)
                {
                    outXC->InsertNextTuple1(inXC->GetComponent(i, 0));
                    outVal->InsertNextTuple1(inXC->GetComponent(i, 1));
                }
            }
            else
            {
                for (int i = 0; i < nPoints; i+=4, index+=4)
                {
                    outXC->InsertNextTuple1( inXC->GetComponent(i,   0));
                    outVal->InsertNextTuple1(inXC->GetComponent(i,   1));
                    outXC->InsertNextTuple1( inXC->GetComponent(i+1, 0));
                    outVal->InsertNextTuple1(inXC->GetComponent(i+1, 1));
                    outXC->InsertNextTuple1( inXC->GetComponent(i+3, 0));
                    outVal->InsertNextTuple1(inXC->GetComponent(i+3, 1));
                    outXC->InsertNextTuple1( inXC->GetComponent(i+2, 0));
                    outVal->InsertNextTuple1(inXC->GetComponent(i+2, 1));
                }
            }
        }
    }

    // 
    // Sort if necessary
    // 
    vtkDataArray *sortedXC;
    vtkDataArray *sortedVal;
    if (requiresSort)
    {
        sortedXC = outXC->NewInstance();
        sortedVal = outVal->NewInstance();
        DoubleIntMap sortedIds;
        nPoints = outXC->GetNumberOfTuples();
        for (int i = 0; i < nPoints; i++)
        {
            x = outXC->GetTuple1(i); 
            sortedIds.insert(DoubleIntMap::value_type(x, i));
        }
        DoubleIntMap::iterator it;
        for (it = sortedIds.begin(); it != sortedIds.end(); it++)
        {
            sortedXC->InsertNextTuple1(outXC->GetTuple1((*it).second));
            sortedVal->InsertNextTuple1(outVal->GetTuple1((*it).second));
        }
        outGrid->SetXCoordinates(sortedXC);
        outGrid->GetPointData()->SetScalars(sortedVal);
        sortedXC->Delete();
        sortedVal->Delete();
    }
    else
    {
        sortedXC  = outXC;
        sortedVal = outVal;
    }
    nPoints = sortedXC->GetNumberOfTuples();
    outGrid->SetDimensions(nPoints, 1, 1);
    outputArray.clear();
    for (int i = 0; i < nPoints; i++)
    {
        outputArray.push_back(sortedXC->GetTuple1(i));
        outputArray.push_back(sortedVal->GetTuple1(i));
    } 

    // Pass a copy of the avtCurveTransform data into the new curve object.
    vtkDataArray *ct = ds[0]->GetFieldData()->GetArray("avtCurveTransform");
    if(ct != 0 && ct->GetNumberOfTuples() == 16)
    {
        ct->Register(NULL);
        outGrid->GetFieldData()->AddArray(ct);
    }

    const char *varname = (pipelineVariable != NULL ? pipelineVariable : "");

    // make sure the outputvar is named.
    sortedVal->SetName(varname);


    //
    //  Clean up.
    // 
    delete [] ds;
    return outGrid;
}

// ****************************************************************************
//  Method: avtCurveConstructorFilter::VerifyInput
//
//  Purpose:
//      Verifies that the input is 2D data, throws an exception if not.
//
//  Programmer: Kathleen Bonnell
//  Creation:   April 26, 2002 
//
// ****************************************************************************
 
void
avtCurveConstructorFilter::VerifyInput(void)
{
    if  (GetInput()->GetInfo().GetAttributes().GetTopologicalDimension() != 1)
    {
        EXCEPTION2(InvalidDimensionsException, "Curve", " Lines ");
    }
}

// ****************************************************************************
//  Method: avtCurveConstructorFilter::ModifyContract
//
//  Purpose:
//    Indicates that we cannot do dynamic load balancing with this filter.  
//
//  Programmer: Kathleen Bonnell
//  Creation:   April 26, 2002 
//
//  Modifications:
//
//    Hank Childs, Tue Feb 19 19:45:43 PST 2008
//    Rename "dynamic" to "streaming", since we really care about whether we
//    are streaming, not about whether we are doing dynamic load balancing.
//    And the two are no longer synonymous.
//
// ****************************************************************************

avtContract_p
avtCurveConstructorFilter::ModifyContract(avtContract_p spec)
{
    spec->NoStreaming();
    return spec;
}

// ****************************************************************************
//  Method: avtCurveConstructorFilter::ModifyContract
//
//  Purpose:
//      Allows the filter to change its output's data object information, which
//      is a description of the data object.
//
//  Programmer: Kathleen Bonnell
//  Creation:   December 23, 2002 
//
//  Modifications:
//    Kathleen Bonnell, Mon Jul 31 16:50:55 PDT 2006
//    Changed Spatial dimension to 1, as curves are not represented as 
//    1D RectilinearGrid.
//
// ****************************************************************************

void
avtCurveConstructorFilter::UpdateDataObjectInfo(void)
{
    GetOutput()->GetInfo().GetAttributes().SetSpatialDimension(1);
}


// ****************************************************************************
//  Method: avtCurveConstructorFilter::PostExecute
//
//  Purpose:
//    Saves the curve data to an output array in PlotInfoAtts.  
//
//  Programmer: Kathleen Bonnell
//  Creation:   June 20, 2006 
//
//  Modifications:
//    Brad Whitlock, Tue Jan  6 16:41:18 PST 2009
//    Store the outputArray into the MapNode data.
//
//    Brad Whitlock, Mon Jul 23 12:15:34 PDT 2012
//    Limit the number of values sent to the client in plot information.
//
// ****************************************************************************

void
avtCurveConstructorFilter::PostExecute(void)
{
#ifdef PARALLEL
    BroadcastDoubleVector(outputArray, PAR_Rank());
#endif
    // Limit outputArray size that we send back to the client.
    if(outputArray.size() < 100000)
    {
        PlotInfoAttributes plotInfoAtts ;
        MapNode data;
        data = outputArray;
        GetOutput()->GetInfo().GetAttributes().AddPlotInformation("Curve", data);
    }
    else
    {
        debug5 << "Curve constructor filter does not send curves that contain "
                  "more than 100K values to the client." << endl;
    }
}

