context: #region Copyright
  /*
  * Copyright (c) 2005-2010, OpenMI Association
  * 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.
  *     * Neither the name of the OpenMI Association 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 "OpenMI Association" ``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 "OpenMI Association" 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.
  */
  #endregion
  #region Copyright
  
  ///////////////////////////////////////////////////////////
  //
  // namespace: org.OpenMI.Utilities.Spatial 
  // purpose: Utility for spatial data-operations within the LinkableComponents wrappers
  // file: ElementMapper.cs
  //
  ///////////////////////////////////////////////////////////
  //
  //    Copyright (C) 2006 OpenMI Association
  //
  //    This library is free software; you can redistribute it and/or
  //    modify it under the terms of the GNU Lesser General Public
  //    License as published by the Free Software Foundation; either
  //    version 2.1 of the License, or (at your option) any later version.
  //
  //    This library is distributed in the hope that it will be useful,
  //    but WITHOUT ANY WARRANTY; without even the implied warranty of
  //    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  //    Lesser General Public License for more details.
  //
  //    You should have received a copy of the GNU Lesser General Public
  //    License along with this library; if not, write to the Free Software
  //    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  //    or look at URL www.gnu.org/licenses/lgpl.html
  //
  //    Contact info: 
  //      URL: www.openmi.org
  //	Email: sourcecode@openmi.org
  //	Discussion forum available at www.sourceforge.net
  //
  //      Coordinator: Roger Moore, CEH Wallingford, Wallingford, Oxon, UK
  //
  ///////////////////////////////////////////////////////////
  //
  //  Original author: Peter Sinding & Jan B. Gregersen, DHI - Water & Environment, Horsholm, Denmark
  //  Created on:      6 April 2005
  //  Version:         1.0.0 
  //
  //  Modification history:
  //  
  //
  ///////////////////////////////////////////////////////////
  
  #endregion
  
  using System;
  using System.Collections;
  using System.Collections.Generic;
  using Oatc.OpenMI.Sdk.Backbone;
  using Oatc.UpwardsComp.Backbone;
  using Oatc.UpwardsComp.Standard;
  using OpenMI.Standard2.TimeSpace;
  using IValueSet=Oatc.UpwardsComp.Standard.IValueSet;
  
  namespace Oatc.UpwardsComp.Spatial
  {
  	/// <summary>
  	/// The ElementMapper converts one ValueSet (inputValues) associated one ElementSet (fromElements)
  	/// to a new ValuesSet (return value of MapValue) that corresponds to another ElementSet 
  	/// (toElements). The conversion is a two step procedure where the first step (Initialize) is 
  	/// executed at initialisation time only, whereas the MapValues is executed during time stepping.
  	/// 
  	/// <p>The Initialize method will create a conversion matrix with the same number of rows as the
  	/// number of elements in the ElementSet associated to the accepting component (i.e. the toElements) 
  	/// and the same number of columns as the number of elements in the ElementSet associated to the 
  	/// providing component (i.e. the fromElements).</p>
  	/// 
  	/// <p>Mapping is possible for any zero-, one- and two-dimensional elemets. Zero dimensional 
  	/// elements will always be points, one-dimensional elements will allways be polylines and two-
  	/// dimensional elements will allways be polygons.</p>
  	/// 
  	/// <p>The ElementMapper contains a number of methods for mapping between the different element types.
  	/// As an example polyline to polygon mapping may be done either as Weighted Mean or as Weighted Sum.
  	/// Typically the method choice will depend on the quantity mapped. Such that state variables such as 
  	/// water level will be mapped using Weighted Mean whereas flux variables such as seepage from river 
  	/// to groundwater will be mapped using Weighted Sum. The list of available methods for a given 
  	/// combination of from and to element types is obtained using the GetAvailableMethods method.</p>
  	/// </summary>
  	public class ElementMapper
  	{
  		private const int NUMBER_OF_AVAILABLE_METHODS = 15;
  		private readonly sMethod[] _availableMethods;
  		private bool _isInitialised;
  		private double[,] _mappingMatrix; // the mapping matrix
  		private int _methodID;
  		private int _numberOfColumns;
  		private int _numberOfRows;
  
  		//=========================================================================================
  		// ElementMapper()
  		//=========================================================================================
  		/// <summary>
  		/// Constructor.
  		/// </summary>
  		public ElementMapper()
  		{
  			_numberOfRows = 0;
  			_numberOfColumns = 0;
  			_isInitialised = false;
  
  			_availableMethods = new sMethod[NUMBER_OF_AVAILABLE_METHODS];
  
  			_availableMethods[0].fromElementsShapeType = ElementType.Point;
  			_availableMethods[0].toElementsShapeType = ElementType.Point;
  			_availableMethods[0].Description = "Nearest";
  			_availableMethods[0].ID = (int) eMethods.PointToPoint.Nearest;
  
  			_availableMethods[1].fromElementsShapeType = ElementType.Point;
  			_availableMethods[1].toElementsShapeType = ElementType.Point;
  			_availableMethods[1].Description = "Inverse";
  			_availableMethods[1].ID = (int) eMethods.PointToPoint.Inverse;
  
  			_availableMethods[2].fromElementsShapeType = ElementType.Point;
  			_availableMethods[2].toElementsShapeType = ElementType.PolyLine;
  			_availableMethods[2].Description = "Nearest";
  			_availableMethods[2].ID = (int) eMethods.PointToPolyline.Nearest;
  
  			_availableMethods[3].fromElementsShapeType = ElementType.Point;
  			_availableMethods[3].toElementsShapeType = ElementType.PolyLine;
  			_availableMethods[3].Description = "Inverse";
  			_availableMethods[3].ID = (int) eMethods.PointToPolyline.Inverse;
  
  			_availableMethods[4].fromElementsShapeType = ElementType.Point;
  			_availableMethods[4].toElementsShapeType = ElementType.Polygon;
  			_availableMethods[4].Description = "Mean";
  			_availableMethods[4].ID = (int) eMethods.PointToPolygon.Mean;
  
  			_availableMethods[5].fromElementsShapeType = ElementType.Point;
  			_availableMethods[5].toElementsShapeType = ElementType.Polygon;
  			_availableMethods[5].Description = "Sum";
  			_availableMethods[5].ID = (int) eMethods.PointToPolygon.Sum;
  
  			_availableMethods[6].fromElementsShapeType = ElementType.PolyLine;
  			_availableMethods[6].toElementsShapeType = ElementType.Point;
  			_availableMethods[6].Description = "Nearest";
  			_availableMethods[6].ID = (int) eMethods.PolylineToPoint.Nearest;
  
  			_availableMethods[7].fromElementsShapeType = ElementType.PolyLine;
  			_availableMethods[7].toElementsShapeType = ElementType.Point;
  			_availableMethods[7].Description = "Inverse";
  			_availableMethods[7].ID = (int) eMethods.PolylineToPoint.Inverse;
  
  			_availableMethods[8].fromElementsShapeType = ElementType.PolyLine;
  			_availableMethods[8].toElementsShapeType = ElementType.Polygon;
  			_availableMethods[8].Description = "Weighted Mean";
  			_availableMethods[8].ID = (int) eMethods.PolylineToPolygon.WeightedMean;
  
  			_availableMethods[9].fromElementsShapeType = ElementType.PolyLine;
  			_availableMethods[9].toElementsShapeType = ElementType.Polygon;
  			_availableMethods[9].Description = "Weighted Sum";
  			_availableMethods[9].ID = (int) eMethods.PolylineToPolygon.WeightedSum;
  
  			_availableMethods[10].fromElementsShapeType = ElementType.Polygon;
  			_availableMethods[10].toElementsShapeType = ElementType.Point;
  			_availableMethods[10].Description = "Value";
  			_availableMethods[10].ID = (int) eMethods.PolygonToPoint.Value;
  
  			_availableMethods[11].fromElementsShapeType = ElementType.Polygon;
  			_availableMethods[11].toElementsShapeType = ElementType.PolyLine;
  			_availableMethods[11].Description = "Weighted Mean";
  			_availableMethods[11].ID = (int) eMethods.PolygonToPolyline.WeightedMean;
  
  			_availableMethods[12].fromElementsShapeType = ElementType.Polygon;
  			_availableMethods[12].toElementsShapeType = ElementType.PolyLine;
  			_availableMethods[12].Description = "Weighted Sum";
  			_availableMethods[12].ID = (int) eMethods.PolygonToPolyline.WeightedSum;
  
  			_availableMethods[13].fromElementsShapeType = ElementType.Polygon;
  			_availableMethods[13].toElementsShapeType = ElementType.Polygon;
  			_availableMethods[13].Description = "Weighted Mean";
  			_availableMethods[13].ID = (int) eMethods.PolygonToPolygon.WeightedMean;
  
  			_availableMethods[14].fromElementsShapeType = ElementType.Polygon;
  			_availableMethods[14].toElementsShapeType = ElementType.Polygon;
  			_availableMethods[14].Description = "Weighted Sum";
  			_availableMethods[14].ID = (int) eMethods.PolygonToPolygon.WeightedSum;
  		}
  
  		//=========================================================================================
  		// Initialise(string methodDescription, IElementSet fromElements,  IElementSet toElements): void
  		//=========================================================================================
  		/// <summary>
  		/// Initialises the ElementMapper. The initialisation includes setting the _isInitialised
  		/// flag and calls UpdateMappingMatrix for claculation of the mapping matrix.
  		/// </summary>
  		///
  		/// <param name="methodDescription">String description of mapping method</param> 
  		/// <param name="fromElements">The IElementSet to map from.</param>
  		/// <param name="toElements">The IElementSet to map to</param>
  		/// 
  		/// <returns>
  		/// The method has no return value.
  		/// </returns>
  		public void Initialise(string methodDescription, IElementSet fromElements, IElementSet toElements)
  		{
  			UpdateMappingMatrix(methodDescription, fromElements, toElements);
  			_isInitialised = true;
  		}
  
  		//=========================================================================================
  		// MapValues(IValueSet inputValues) : IValueSet
  		//=========================================================================================
  		/// <summary>
  		/// MapValues calculates a IValueSet through multiplication of an inputValues IValueSet
  		/// vector or matrix (ScalarSet or VectorSet) on to the mapping maprix. IScalarSets maps
  		/// to IScalarSets and IVectorSets maps to IVectorSets.
  		/// </summary>
  		/// 
  		/// <remarks>
  		/// Mapvalues is called every time a georeferenced link is evaluated.
  		/// </remarks>
  		/// 
  		/// <param name="inputValues">IValueSet of values to be mapped.</param>
  		/// 
  		/// <returns>
  		/// A IValueSet found by mapping of the inputValues on to the toElementSet.
  		/// </returns>
  		public IValueSet MapValues(IValueSet inputValues)
  		{
  			if (!_isInitialised)
  			{
  				throw new Exception("ElementMapper objects needs to be initialised before the MapValue method can be used");
  			}
  			if (!inputValues.Count.Equals(_numberOfColumns))
  			{
  				throw new Exception("Dimension mismatch between inputValues and mapping matrix");
  			}
  			if (inputValues is IScalarSet)
  			{
  				double[] outValues = new double[_numberOfRows];
  				//--- Multiply the Values vector with the MappingMatrix ---
  				for (int i = 0; i < _numberOfRows; i++)
  				{
  					outValues[i] = 0;
  					for (int n = 0; n < _numberOfColumns; n++)
  					{
  						outValues[i] += _mappingMatrix[i, n]*((IScalarSet) inputValues).GetScalar(n); //(remove)inValues[n];
  					}
  				}
  				ScalarSet outputValues = new ScalarSet(outValues);
  				return outputValues;
  			}
  /*
  			else if (inputValues is IVectorSet)
  			{
  				Vector[] outValues = new Vector[_numberOfRows];
  				//--- Multiply the Values vector with the MappingMatrix ---
  				for (int i = 0; i < _numberOfRows; i++)
  				{
  					outValues[i].XComponent = 0;
  					outValues[i].YComponent = 0;
  					outValues[i].ZComponent = 0;
  					for (int n = 0; n < _numberOfColumns; n++)
  					{
  						outValues[i].XComponent += _mappingMatrix[i, n]*((IVectorSet) inputValues).GetVector(n).XComponent;
  						outValues[i].YComponent += _mappingMatrix[i, n]*((IVectorSet) inputValues).GetVector(n).YComponent;
  						outValues[i].ZComponent += _mappingMatrix[i, n]*((IVectorSet) inputValues).GetVector(n).ZComponent;
  					}
  				}
  				VectorSet outputValues = new VectorSet(outValues);
  				return outputValues;
  			}
  */
  			else
  			{
  				throw new Exception("Invalid datatype used for inputValues parameter. MapValues failed");
  			}
  		}
  
  		//=========================================================================================
  		//UpdateMappingMatrix(string methodDescription, IElementSet fromElements, IElementSet toElements) : void
  		//=========================================================================================
  		/// <summary>
  		/// Calculates the mapping matrix between fromElements and toElements. The mapping method 
  		/// is decided from the combination of methodDescription, fromElements.ElementType and 
  		/// toElements.ElementType. 
  		/// The valid values for methodDescription is obtained through use of the 
  		/// GetAvailableMethods method.
  		/// </summary>
  		/// 
  		/// <remarks>
  		/// UpdateMappingMatrix is called during initialisation. UpdateMappingMatrix must be called prior
  		/// to Mapvalues.
  		/// </remarks>
  		/// 
  		/// <param name="methodDescription">String description of mapping method</param> 
  		/// <param name="fromElements">The IElementset to map from.</param>
  		/// <param name="toElements">The IElementset to map to</param>
  		///
  		/// <returns>
  		/// The method has no return value.
  		/// </returns>
  		public void UpdateMappingMatrix(string methodDescription, IElementSet fromElements, IElementSet toElements)
  		{
  			try
  			{
  				ElementSetChecker.CheckElementSet(fromElements);
  				ElementSetChecker.CheckElementSet(toElements);
  
  				_methodID = GetMethodID(methodDescription, fromElements.ElementType, toElements.ElementType);
  				_numberOfRows = toElements.ElementCount;
  				_numberOfColumns = fromElements.ElementCount;
  				_mappingMatrix = new double[_numberOfRows,_numberOfColumns];
  
  				if (fromElements.ElementType == ElementType.Point && toElements.ElementType == ElementType.Point)
  					// Point to Point
  				{
  					#region
  
  					try
  					{
  						for (int i = 0; i < _numberOfRows; i++)
  						{
  							XYPoint ToPoint = CreateXYPoint(toElements, i);
  							for (int j = 0; j < _numberOfColumns; j++)
  							{
  								XYPoint FromPoint = CreateXYPoint(fromElements, j);
  								_mappingMatrix[i, j] = XYGeometryTools.CalculatePointToPointDistance(ToPoint, FromPoint);
  							}
  						}
  
  						if (_methodID.Equals((int) eMethods.PointToPoint.Nearest))
  						{
  							for (int i = 0; i < _numberOfRows; i++)
  							{
  								double MinDist = _mappingMatrix[i, 0];
  								for (int j = 1; j < _numberOfColumns; j++)
  								{
  									if (_mappingMatrix[i, j] < MinDist)
  									{
  										MinDist = _mappingMatrix[i, j];
  									}
  								}
  								int Denominator = 0;
  								for (int j = 0; j < _numberOfColumns; j++)
  								{
  									if (_mappingMatrix[i, j] == MinDist)
  									{
  										_mappingMatrix[i, j] = 1;
  										Denominator++;
  									}
  									else
  									{
  										_mappingMatrix[i, j] = 0;
  									}
  								}
  								for (int j = 0; j < _numberOfColumns; j++)
  								{
  									_mappingMatrix[i, j] = _mappingMatrix[i, j]/Denominator;
  								}
  							}
  						} // if (_methodID.Equals((int) eMethods.PointToPoint.Nearest))   
  						else if (_methodID.Equals((int) eMethods.PointToPoint.Inverse))
  						{
  							for (int i = 0; i < _numberOfRows; i++)
  							{
  								double MinDist = _mappingMatrix[i, 0];
  								for (int j = 1; j < _numberOfColumns; j++)
  								{
  									if (_mappingMatrix[i, j] < MinDist)
  									{
  										MinDist = _mappingMatrix[i, j];
  									}
  								}
  								if (MinDist == 0)
  								{
  									int Denominator = 0;
  									for (int j = 0; j < _numberOfColumns; j++)
  									{
  										if (_mappingMatrix[i, j] == MinDist)
  										{
  											_mappingMatrix[i, j] = 1;
  											Denominator++;
  										}
  										else
  										{
  											_mappingMatrix[i, j] = 0;
  										}
  									}
  									for (int j = 0; j < _numberOfColumns; j++)
  									{
  										_mappingMatrix[i, j] = _mappingMatrix[i, j]/Denominator;
  									}
  								}
  								else
  								{
  									double Denominator = 0;
  									for (int j = 0; j < _numberOfColumns; j++)
  									{
  										_mappingMatrix[i, j] = 1/_mappingMatrix[i, j];
  										Denominator = Denominator + _mappingMatrix[i, j];
  									}
  									for (int j = 0; j < _numberOfColumns; j++)
  									{
  										_mappingMatrix[i, j] = _mappingMatrix[i, j]/Denominator;
  									}
  								}
  							}
  						} // else if (_methodID.Equals((int) eMethods.PointToPoint.Inverse))
  						else
  						{
  							throw new Exception("methodDescription unknown for point point mapping");
  						}
  							// else if (_methodID.Equals((int) eMethods.PointToPoint.Nearest)) and else if (_methodID.Equals((int) eMethods.PointToPoint.Inverse))
  					}
  					catch (Exception e) // Catch for all of the Point to Point part
  					{
  						throw new Exception("Point to point mapping failed", e);
  					}
  
  					#endregion
  				}
  				else if (fromElements.ElementType == ElementType.Point && toElements.ElementType == ElementType.PolyLine)
  					// Point to PolyLine
  				{
  					#region
  
  					try
  					{
  						for (int i = 0; i < _numberOfRows; i++)
  						{
  							XYPolyline toPolyLine = CreateXYPolyline(toElements, i);
  							for (int j = 0; j < _numberOfColumns; j++)
  							{
  								XYPoint fromPoint = CreateXYPoint(fromElements, j);
  								_mappingMatrix[i, j] = XYGeometryTools.CalculatePolylineToPointDistance(toPolyLine, fromPoint);
  							}
  						}
  
  						if (_methodID.Equals((int) eMethods.PointToPolyline.Nearest))
  						{
  							for (int i = 0; i < _numberOfRows; i++)
  							{
  								double MinDist = _mappingMatrix[i, 0];
  								for (int j = 1; j < _numberOfColumns; j++)
  								{
  									if (_mappingMatrix[i, j] < MinDist)
  									{
  										MinDist = _mappingMatrix[i, j];
  									}
  								}
  								int denominator = 0;
  								for (int j = 0; j < _numberOfColumns; j++)
  								{
  									if (_mappingMatrix[i, j] == MinDist)
  									{
  										_mappingMatrix[i, j] = 1;
  										denominator++;
  									}
  									else
  									{
  										_mappingMatrix[i, j] = 0;
  									}
  								}
  								for (int j = 0; j < _numberOfColumns; j++)
  								{
  									_mappingMatrix[i, j] = _mappingMatrix[i, j]/denominator;
  								}
  							}
  						} // if (_methodID.Equals((int) eMethods.PointToPolyline.Nearest))
  						else if (_methodID.Equals((int) eMethods.PointToPolyline.Inverse))
  						{
  							for (int i = 0; i < _numberOfRows; i++)
  							{
  								double minDist = _mappingMatrix[i, 0];
  								for (int j = 1; j < _numberOfColumns; j++)
  								{
  									if (_mappingMatrix[i, j] < minDist)
  									{
  										minDist = _mappingMatrix[i, j];
  									}
  								}
  								if (minDist == 0)
  								{
  									int denominator = 0;
  									for (int j = 0; j < _numberOfColumns; j++)
  									{
  										if (_mappingMatrix[i, j] == minDist)
  										{
  											_mappingMatrix[i, j] = 1;
  											denominator++;
  										}
  										else
  										{
  											_mappingMatrix[i, j] = 0;
  										}
  									}
  									for (int j = 0; j < _numberOfColumns; j++)
  									{
  										_mappingMatrix[i, j] = _mappingMatrix[i, j]/denominator;
  									}
  								}
  								else
  								{
  									double denominator = 0;
  									for (int j = 0; j < _numberOfColumns; j++)
  									{
  										_mappingMatrix[i, j] = 1/_mappingMatrix[i, j];
  										denominator = denominator + _mappingMatrix[i, j];
  									}
  									for (int j = 0; j < _numberOfColumns; j++)
  									{
  										_mappingMatrix[i, j] = _mappingMatrix[i, j]/denominator;
  									}
  								}
  							}
  						} // else if (_methodID.Equals((int) eMethods.PointToPolyline.Inverse))
  						else // if _methodID != Nearest and Inverse
  						{
  							throw new Exception("methodDescription unknown for point to polyline mapping");
  						}
  					}
  					catch (Exception e) // Catch for all of the Point to Polyline part 
  					{
  						throw new Exception("Point to polyline mapping failed", e);
  					}
  
  					#endregion
  				}
  				else if (fromElements.ElementType == ElementType.Point && toElements.ElementType == ElementType.Polygon)
  					// Point to Polygon
  				{
  					#region
  
  					try
  					{
  						for (int i = 0; i < _numberOfRows; i++)
  						{
  							XYPolygon polygon = CreateXYPolygon(toElements, i);
  							int count = 0;
  							XYPoint point;
  							for (int n = 0; n < _numberOfColumns; n++)
  							{
  
gt: 								point = CreateXYPoint(fromElements, n);