Changeset 58f6e15 in tol for tol_pkg


Ignore:
Timestamp:
Nov 12, 2009, 9:36:20 PM (17 years ago)
Author:
Víctor de Buen <vdebuen@…>
Branches:
1941-docker-image-to-test, devel, master
Children:
8f0fb94
Parents:
ae5159f
Message:

Moving matrix methods to new package MatQuery and creating GIS related classes

git-svn-id: ​https://www.tol-project.org/svn/tolp/trunk@1696 be4846ae-f051-483a-831f-b64e7d08527f

File:
1 edited

Legend:

Unmodified
Added
Removed

  • tol_pkg/GisTools/GisTools.tol

    rae5159f r58f6e15  
    22//////////////////////////////////////////////////////////////////////////////
    33// FILE   : GisTools.tol
    4 // PURPOSE: Defines NameBlock GisTools
    5 //////////////////////////////////////////////////////////////////////////////
    6 
     4// PURPOSE: Defines Package GisTools
     5//////////////////////////////////////////////////////////////////////////////
     6
     7//////////////////////////////////////////////////////////////////////////////
    78NameBlock GisTools = [[
     9//////////////////////////////////////////////////////////////////////////////
     10
     11#Require MatQuery;
    812
    913Text _.autodoc.description = "Functions and classes related with "
    10 "Geographic Information System (GIS), selection and classification methods "
    11 "over matricial data over massive geographical points.";
     14"Geographic Information System (GIS), with all type of methods "
     15"over large matricial geographical points.";
    1216Text _.autodoc.url = "";
    13 Set _.autodoc.keys = [["GIS","Geographic","select","classify"]];
     17Set _.autodoc.keys = [["GIS","Geographic","point","distance"]];
    1418Set _.autodoc.authors = [[
    1519   "vdebuen@tol-project.org"   ]];
    1620
    1721//////////////////////////////////////////////////////////////////////////////
    18 Text _.autodoc.member.SelectRowsWithValue =
    19 "Selects rows of a column matrix whose cell is exactly the specified value";
    20 Set SelectRowsWithValue(Matrix M, Real value)
    21 //////////////////////////////////////////////////////////////////////////////
    22 {
    23   Real r = Rows(M);
    24   VMatrix V = Rand(r,1,value,value);
    25   VMatrix vM = Mat2VMat(M);
    26   VMatrix match = EQ(vM,V);
    27   Matrix select = VMat2Triplet(Convert(match,"Cholmod.R.Triplet"));
    28   If(!Rows(select),Copy(Empty), MatSet(Tra(SubCol(select,[[1]])))[1])
    29 };
    30 
    31 
    32 //////////////////////////////////////////////////////////////////////////////
    33 Text _.autodoc.member.SelectRowsInInterval =
    34 "Selects rows of a column matrix whose cell is inside the interval defined "
    35 "as: \n"
    36 "  [min,max]: if minIncluded=True  and minIncluded=True\n"
    37 "  (min,max]: if minIncluded=False and minIncluded=True\n"
    38 "  [min,max): if minIncluded=True  and minIncluded=False\n"
    39 "  (min,max): if minIncluded=False and minIncluded=False\n";
    40 Set SelectRowsInInterval(
    41  Matrix M,
    42  Real minIncluded,
    43  Real min,
    44  Real max,
    45  Real maxIncluded )
    46 //////////////////////////////////////////////////////////////////////////////
    47 {
    48   Real r = Rows(M);
    49   VMatrix Min = Rand(r,1,min,min);
    50   VMatrix Max = Rand(r,1,max,max);
    51   VMatrix vM = Mat2VMat(M);
    52   VMatrix match =
    53   Case(
    54      minIncluded &  maxIncluded, And(LE(Min,vM),LE(vM,Max)),
    55     !minIncluded &  maxIncluded, And(LT(Min,vM),LE(vM,Max)),
    56      minIncluded & !maxIncluded, And(LE(Min,vM),LT(vM,Max)),
    57     !minIncluded & !maxIncluded, And(LT(Min,vM),LT(vM,Max))
    58   );
    59   Matrix select = VMat2Triplet(Convert(match,"Cholmod.R.Triplet"));
    60   If(!Rows(select),Copy(Empty), MatSet(Tra(SubCol(select,[[1]])))[1])
    61 };
    62 
    63 //////////////////////////////////////////////////////////////////////////////
    64 Text _.autodoc.member.SelectRowsInHyperCube = "Selects rows of matrix "
    65 "M that are inside the hipercube whose side length is 2*d and is centered in "
    66 "point given by row matrix P ";
    67 Set SelectRowsInHyperCube(Matrix M, Matrix P, Real d)
    68 //////////////////////////////////////////////////////////////////////////////
    69 {
    70   Real r = Rows(M);
    71   Real c = Columns(M);
    72   Matrix D = Rand(r,1,d,d);
    73   VMatrix match = Group("And",For(1,c, VMatrix(Real k)
    74   {
    75     Mat2VMat(LE(Abs(SubCol(M,[[k]])-MatDat(P,1,k)),D))
    76   }));
    77   Matrix select = VMat2Triplet(Convert(match,"Cholmod.R.Triplet"));
    78   If(!Rows(select),Copy(Empty), MatSet(Tra(SubCol(select,[[1]])))[1])
    79 };
    80 
    81 //////////////////////////////////////////////////////////////////////////////
    82 Text _.autodoc.member.SelectRowsInHyperSphere = "Selects rows of "
    83 "matrix M that are inside the hipersphere with specified  ratio and is "
    84 "centered in point given by row matrix P";
    85 Set SelectRowsInHyperSphere(Matrix M, Matrix P, Real ratio)
    86 //////////////////////////////////////////////////////////////////////////////
    87 {
    88   Real r = Rows(M);
    89   Real c = Columns(M);
    90   Matrix R = Rand(r,1,ratio,ratio);
    91   Matrix dist = Sqrt(SetSum(For(1,c, Matrix(Real k)
    92   {
    93     RPow(SubCol(M,[[k]])-MatDat(P,1,k),2)
    94   })));
    95   VMatrix match = Mat2VMat(LE(dist,R));
    96   Matrix select = VMat2Triplet(Convert(match,"Cholmod.R.Triplet"));
    97   If(!Rows(select),Copy(Empty), MatSet(Tra(SubCol(select,[[1]])))[1])
    98 };
    99 
    100 //////////////////////////////////////////////////////////////////////////////
    101 Text _.autodoc.member.ClassifyRowsByRank = "Classifies the rows of a column "
    102 "matrix by the ranking position respect a set of R sorted numbers. The class "
    103 "of the i-th row will be assigned in this way:\n"
    104 " class[i] = 1   <=>                 M[i] < rankBound[1]\n"
    105 " class(i] = k+1 <=> rankBound[k] <= M[i] < rankBound[k+1]; k=1...R\n"
    106 " class[i] = R+1 <=> rankBound[R] <= M[i]\n";
    107 Matrix ClassifyRowsByRank(Matrix M, Set rankBound)
    108 //////////////////////////////////////////////////////////////////////////////
    109 {
    110   Real m = Rows(M);
    111   Real R = Card(rankBound);
    112   Real min =
    113   {
    114     Real aux = MatMin(M);
    115     aux - Abs(aux) - 1
    116   };
    117   Real max =
    118   {
    119     Real aux = MatMax(M);
    120     aux + Abs(aux) + 1
    121   };
    122   Set rnk = SetOfReal(min)<<rankBound<<SetOfReal(max);
    123   VMatrix class = Zeros(m,1);
    124   VMatrix vM = Mat2VMat(M);
    125   Real r = 1;
    126   VMatrix constant = Rand(m,1,1,1);
    127   While(r<=R+1,
    128   {
    129     VMatrix left.match = LE(constant*rnk[r],vM);
    130     VMatrix right.match = LT(vM,constant*rnk[r+1]);
    131     VMatrix class := class + And(left.match,right.match)*r;
    132     Real r := r+1
    133   });
    134   VMat2Mat(class)
     22//Abstract class to handle with a bidimensional point.
     23Class @Point2D
     24//////////////////////////////////////////////////////////////////////////////
     25{
     26  ////////////////////////////////////////////////////////////////////////////
     27  Static Text _.autodoc.member.numberOfPoints =
     28  "Returns the number of stored points";
     29  Real numberOfPoints(Real void);
     30  ////////////////////////////////////////////////////////////////////////////
     31
     32  ////////////////////////////////////////////////////////////////////////////
     33  Static Text _.autodoc.member.x =
     34  "Returns the first coordinate of the k-th point";
     35  Real x(Real k);
     36  ////////////////////////////////////////////////////////////////////////////
     37
     38  ////////////////////////////////////////////////////////////////////////////
     39  Static Text _.autodoc.member.y =
     40  "Returns the second coordinate of the k-th point";
     41  Real y(Real k)
     42  ////////////////////////////////////////////////////////////////////////////
     43};
     44
     45//////////////////////////////////////////////////////////////////////////////
     46//Class to handle with just one bidimensional point that
     47//isnŽt externally stored
     48Class @Point2D.Simple : @Point2D
     49//////////////////////////////////////////////////////////////////////////////
     50{
     51  Real _.x;
     52  Real _.y;
     53
     54  //Just one point is stored
     55  Real numberOfPoints(Real void) { 1 };
     56  //Argument k is unused
     57  Real x(Real k) { _.x };
     58  //Argument k is unused
     59  Real y(Real k) { _.y };
     60  ////////////////////////////////////////////////////////////////////////////
     61  Static Text _.autodoc.member.New = "Creates an instance of @Point2D.Simple";
     62  @Point2D.Simple New(Real x, Real y)
     63  ////////////////////////////////////////////////////////////////////////////
     64  {
     65    @Point2D.Simple new = [[
     66      Real _.x = x;
     67      Real _.y = y ]];
     68    new
     69  }
     70};
     71
     72//////////////////////////////////////////////////////////////////////////////
     73//Class to store al large set of points as rows of a matrix with 2 columns.
     74Class @Point2D.Stored
     75//////////////////////////////////////////////////////////////////////////////
     76{
     77  VMatrix _.xy;
     78
     79  Real numberOfPoints(Real void) { VRows(_.xy) };
     80  Real x(Real k) { VMatDat(_.xy,k,1) };
     81  Real y(Real k) { VMatDat(_.xy,k,1) };
     82
     83  @Point2D.Simple center(Real void)
     84  {
     85    Real m = VRows(_.xy);
     86    Matrix c = Rand(m,1,1/m,1/m)*_.xy;
     87    @Point2D.Simple::New(MatDat(c,1,1), MatDat(c,1,2))
     88  };
     89
     90  ////////////////////////////////////////////////////////////////////////////
     91  Static Text _.autodoc.member.New = "Creates an instance of @Point2D.Stored";
     92  @Point2D.Stored New(VMatrix xy)
     93  ////////////////////////////////////////////////////////////////////////////
     94  {
     95    If(Colunms(xy)>2,
     96      WriteLn("[@Point2D.Stored::New] Only 2 first of "<<Colunms(xy)+
     97      " columns of argument xy will be used","W"));
     98    If(Colunms(xy)<2,
     99      WriteLn("[@Point2D.Stored::New] Almost 2 columns instead of "<<
     100      Colunms(xy)+" are needed for argument xy","E"));
     101    @Point2D.Simple new = [[
     102      VMatrix _.xy = xy ]];
     103    new
     104  }
     105};
     106
     107//////////////////////////////////////////////////////////////////////////////
     108//Abstract class to inherite bidimensional interpolating classes.
     109Class @Interpolate2D
     110//////////////////////////////////////////////////////////////////////////////
     111{
     112  //Bidimensional points in cartesian coordinates. Dim: _.m x 2
     113  @Point2D.Stored _.P;
     114  //Tarjet data to be interpolated, Dim : _.m x _.n
     115  VMatrix _.F;
     116
     117  ////////////////////////////////////////////////////////////////////////////
     118  //Auxiliar members
     119  ////////////////////////////////////////////////////////////////////////////
     120  //Number of interpolation points
     121  Real _.m = ?;
     122  //Number of independent interpolation functions.
     123  //When _.n==1 itŽs called the scalar interpolating case
     124  //When _.n>=2 itŽs called the vectorial interpolating case
     125  Real _.n = ?;
     126
     127  ////////////////////////////////////////////////////////////////////////////
     128  Static Text _.autodoc.member.build =
     129  "Initializes all members which are needed to build the interpolating "
     130  "methods";
     131  Real build(Real void);
     132  ////////////////////////////////////////////////////////////////////////////
     133
     134  ////////////////////////////////////////////////////////////////////////////
     135  Static Text _.autodoc.member.eval.vector =
     136  "Evaluates all interpolation function at point x,y";
     137  VMatrix eval.vector(@Point2D p);
     138  ////////////////////////////////////////////////////////////////////////////
     139
     140  ////////////////////////////////////////////////////////////////////////////
     141  Static Text _.autodoc.member.eval =
     142  "Evaluates the first interpolation function at point x,y. It's useful in "
     143  "the scalar case, when _.F has just one column (_.n==1)";
     144  Real eval(@Point2D p)
     145  ////////////////////////////////////////////////////////////////////////////
     146  {
     147    VMatDat(eval.vector(p),1,1)
     148  }
     149};
     150
     151//////////////////////////////////////////////////////////////////////////////
     152//Abstract class to inherite bidimensional weighted interpolating classes.
     153Class @Interpolate2D.Weighted : @Interpolate2D
     154//////////////////////////////////////////////////////////////////////////////
     155{
     156  //In order to be more efficient, the points farther of this distance will
     157  //be skipped due its weight is almost null
     158  Real _.maxDistance;
     159
     160  ////////////////////////////////////////////////////////////////////////////
     161  //Auxiliar members
     162  ////////////////////////////////////////////////////////////////////////////
     163  Set     _.x.range   = Copy(Empty);
     164  VMatrix _.x.class   = Rand(0,0,0,0);
     165  Set     _.y.range   = Copy(Empty);
     166  VMatrix _.y.class   = Rand(0,0,0,0);
     167  Set     _.point.idx = Rand(0,0,0,0);
     168  Set     _.grid      = Copy(Empty);
     169 
     170  ////////////////////////////////////////////////////////////////////////////
     171  Real build(Real void)
     172  ////////////////////////////////////////////////////////////////////////////
     173  {
     174    Real _.m := VRows(_.P::_.XY);
     175    Real _.n := VRows(_.F);
     176    If(VColunms(_.P::_.XY)!=2,
     177      WriteLn("Interpolate2D.Weighted::_.P::_.XY must have 2 columns instead "
     178              "of "<<Colunms(_.P::_.XY)));
     179    If(VRows(_.F)!=m,
     180      WriteLn("Interpolate2D.Weighted::_.XY must have 2 columns instead of "
     181              <<Colunms(_.P::_.XY)));
     182
     183    Set _.point.idx := Range(1,_.m,1);
     184    VMatrix x = SubCol(_.P::_.xy,[[1]]);
     185    VMatrix y = SubCol(_.P::_.xy,[[2]]);
     186    Set _.x.range := Range(VMatMin(x),VMatMax(x),_.maxDistance);
     187    Set _.y.range := Range(VMatMin(y),VMatMax(y),_.maxDistance);
     188    VMatrix _.x.class := MatQuery::ClassifyRowsByRank(x,_.x.range);
     189    VMatrix _.y.class := MatQuery::ClassifyRowsByRank(y,_.y.range);
     190    Set xy.class = Classify(point.idx, Real(Real a, Real b)
     191    {
     192      Real x.cmp = Compare(VMatDat(_.x.class,a,1),VMatDat(_.x.class,b,1));
     193      If(x.cmp, x.cmp,
     194      Real y.cmp = Compare(VMatDat(_.y.class,a,1),VMatDat(_.y.class,b,1)))
     195    });
     196    Set _.grid := For(1,Card(_.x.range)+1,Set(Real i)
     197    {
     198      For(1,Card(_.y.range)+1,Set(Real j)
     199      {
     200        Copy(Empty)
     201      })
     202    });
     203    Set EvalSet(xy.class, Real(Set cl)
     204    {
     205      Real k = cl[1];
     206      Real x.class = MatDat(_.x.class,k,1);
     207      Real y.class = MatDat(_.y.class,k,1);
     208      Set _.grid[x.class][y.class] := cl;
     209      k
     210    });
     211    True
     212  };
     213
     214  ////////////////////////////////////////////////////////////////////////////
     215  Static Text _.autodoc.member.get.distance =
     216  "Evaluates the weight between a given point p and each point of a 2 "
     217  "columns matrix of coordinates";
     218  VMatrix get.distance(@Point2D p, VMatrix xy)
     219  ////////////////////////////////////////////////////////////////////////////
     220  {
     221    Real m = VRows(xy);
     222    VMatrix p.x = Rand(m,1,p::X(1),p::X(1));
     223    VMatrix p.y = Rand(m,1,p::Y(1),p::Y(1));
     224    VMatrix x = SubCol(xy,[[1]]);
     225    VMatrix y = SubCol(xy,[[1]]);
     226    RPow(RPow(x,2)+RPow(y,2),0.5)
     227  };
     228 
     229  ////////////////////////////////////////////////////////////////////////////
     230  VMatrix eval.vector(@Point2D p)
     231  ////////////////////////////////////////////////////////////////////////////
     232  {
     233    VMatrix p.x = Rand(1,1,p::X(1),p::X(1));
     234    VMatrix p.y = Rand(1,1,p::Y(1),p::Y(1));
     235    VMatrix p.xy = p.x | p.y;
     236    Real p.x.class = VMatDat(MatQuery::ClassifyRowsByRank(p.x,_.x.range),1,1);
     237    Real p.y.class = VMatDat(MatQuery::ClassifyRowsByRank(p.y,_.y.range),1,1);
     238    Real i.min = Max(1,p.x.class-1);
     239    Real i.max = Min(Card(_.x.range)+1,p.x.class+1);
     240    Real j.min = Max(1,p.y.class-1);
     241    Real j.max = Min(Card(_.y.range)+1,p.y.class+1);
     242    Set p.candidates = Bingroup("<<",For(i.min, i.max, Real(Real i)
     243    {
     244      Bingroup("<<",For(j.min, j.max, Real(Real j)
     245      {
     246        _.grid[i][j] 
     247      }))
     248    }));
     249    VMatrix xy.candidates = SubRow(_.P::_.xy, p.neighborhood);
     250    VMatrix d.candidates = get.distance(p,xy.candidates);
     251    Set neighborhood = MatQuery::SelectRowsInInterval(d.candidates,
     252      False, 0, _.maxDistance, True);
     253    VMatrix d = SubRows(d.candidates, neighborhood);
     254    VMatrix w = RPow(d,-0.5);
     255    Tra(w/VMatSum(w)) * SubRow(SubRow(_.F,p.candidates),neighborhood)
     256  };
     257
     258  Static Text _.autodoc.member.New = "Creates an instance of "
     259  "@Interpolate2D.Weighted";
     260  @Interpolate2D.Weighted New(VMatrix XY, VMatrix F, Real maxDistance)
     261  {
     262    @Interpolate2D.Weighted new =
     263    [[
     264      @Point2D.Stored _.P = @Point2D.Stored::New(XY);
     265      VMatrix _.F = F;
     266      Real _.maxDistance = maxDistance
     267    ]]
     268  }
    135269}
    136270
    137 
    138271]];
     272
     273
Note: See TracChangeset for help on using the changeset viewer.