close Warning: Can't synchronize with repository "(default)" (/var/svn/mms does not appear to be a Subversion repository.). Look in the Trac log for more information.

Changes between Version 3 and Version 5 of Ticket #79


Ignore:
Timestamp:
Aug 25, 2010, 7:22:12 AM (15 years ago)
Author:
Pedro Gea
Comment:

Legend:

Unmodified
Added
Removed
Modified

  • Ticket #79

    • Property Priority changed from critical to blocker
    • Property Status changed from new to accepted
    • Property Owner changed from josp to Pedro Gea
    • Property Component changed from General to Tests

  • Ticket #79 – Description

    v3 v5  
    1 Comprobar si existe alguna combinación lineal entre parámetros de un nodo
    2 observacional, teniendo en cuenta si éste tiene o no una distribución a
    3 priori y/o está dentro de una jerarquía. Y además debe informar cuál es la
    4 combinación lineal.
     1Comprobar si existe alguna combinación lineal entre parámetros de un nodo observacional, teniendo en cuenta si éste tiene o no una distribución a priori y/o está dentro de una jerarquía. Y además debe informar cuál es la combinación lineal.
    52
    6 Os escribo debajo las siete funciones que hacen este proceso en bsrlayer,
    7 una llama a las otras, pero usan sets de información que no existen en
    8 mms, si decidís usarlas como guía preguntadme cualquier duda que tengáis.
    9 Al final, si detecta una combinación lineal, escribe los parámetros que la
    10 forman en un fichero log, en forma de query.
     3Se adjuntan las siete funciones que hacen este proceso en BSRLayer, una llama a las otras, pero usan sets de información que no existen en MMS, si decidís usarlas como guía preguntadme cualquier duda que tengáis.
     4Al final, si detecta una combinación lineal, escribe los parámetros que la forman en un fichero log, en forma de query.
    115
    12 Una cosa a tener en cuenta. El proceso en general funcionaba en torno a la
    13 función SVD de tol, sin embargo, como esta función no es exacta sino
    14 aproximada, se probaron los diferentes métodos: Golub_Reinsch, Jacobi y
    15 Sparse. De todos ellos, el más fiable fué el primero. El Sparse era muy
    16 rápido cuando las matrices a factorizar tenían densidades pequeñas (<
    17 0.01), y el Jacobi mostraba una salto de nivel muy claro entre los
    18 autovalores de las variables que no son combinaciones lineales de los que
    19 sí lo son. Sin embargo éstos dos últimos fallaban a veces, detectando
    20 combinaciones lineales donde no las había o al revés. Por defecto se
    21 decidió usar Golub_Reinsch con un redondeo en 10-10
    22 
    23 {{{
    24 //////////////////////////////////////////////////////////////////////////////
    25 Real CheckLinCom( Set OutputInfo, Set InputsInfo, Set AllHierarchyInfo,
    26   Text obs.node, Real vl_inf, Real vl_pri, Text logRoute)
    27 //////////////////////////////////////////////////////////////////////////////
    28 {
    29   // We extract all parameters involved in any hierarchy for node obs.node
    30   Set sHieNodeInfo = BinGroup("<<",EvalSet(AllHierarchyInfo,Set (Set s)
    31   {
    32     Set prev1 = Extract(s,2)|Extract(s,1);
    33 
    34     Set prev2 = Select(prev1,Real (Set x)
    35     {
    36       Real If(x[1] == obs.node, True, False)
    37     }) - [[ Empty ]];
    38 
    39     Set BinGroup("<<",Extract(prev2,2))
    40   }));
    41 
    42 
    43   Real c = 0;
    44 
    45   Set inputs = Select(OutputInfo::InputDB,Real (Anything x)
    46   {
    47     Real c := c+1;
    48 
    49     // Real flag is a label that tell us if an input has not any prior, or
    50     // has not any unknown value, or if it is involved in any hierarchy
    51 tree
    52 
    53     Real flag =
    54       If(
    55         Name(InputsInfo[c]) <: sHieNodeInfo,
    56         Real False,
    57         If(
    58           Not(And( BinEQ( ?, (InputsInfo[c])->prior_mu ),
    59                    BinEQ( ?, (InputsInfo[c])->prior_sigma ))),
    60           Real False,
    61           If(
    62             Grammar(x) == "Serie",
    63             Real Not(MatSum(IsUnknown(SerMat(x)))),// HasUnknown is slower
    64             If(
    65               Grammar(x) == "Matrix",
    66               Real Not(MatSum(IsUnknown(x))),// HasUnknown is slower
    67               Real True
    68             )
    69           )
    70         )
    71       );
    72 
    73     flag
    74   });
    75 
    76   Real If( Not(Card(inputs)),
    77     Real True,
    78     Real LinComWarning(OutputInfo,inputs,logRoute,id_model,obs.node)
    79   )
    80 };
    81 //////////////////////////////////////////////////////////////////////////////
    82 PutDescription(
    83 "It checks a set of input series (or matrix) to find linear combinations
    84 between the observational node parameters. It needs an OutputInfo BSR set
    85 of
    86 input series, the observational node, the real variable vl_pri (True or
    87 False)
    88 and a Path to store a log file with the report of the parameters that has
    89 linear combinations between them.
    90 If exist linear combination, it writes in log file and kill TOL",
    91 CheckLinCom);
    92 //////////////////////////////////////////////////////////////////////////////
    93 
    94 
    95 //////////////////////////////////////////////////////////////////////////////
    96 Real LinComWarning(Set OutputInfo, Set inputs, Text logRoute, Text
    97 id_model,
    98   Text node)
    99 //////////////////////////////////////////////////////////////////////////////
    100 {
    101   Real matrix = If( id_dating == "Matrix",True,False);
    102 
    103   Set lcw = If( Not(matrix), LinComWarningSerie( OutputInfo,inputs ),
    104     LinComWarningMatrix( inputs ));
    105 
    106   Real If( EQ(Card(lcw),0), True,
    107   {
    108     Text WriteLn("[LinComWarning] THERE ARE "+IntText(Card(lcw))+" LINEAR
    109 "+
    110       "COMBINATIONS OF INPUT SERIES CHECK LOG DIRECTORY!!", Text "E");
    111 
    112     Text txt2File =
    113     "update "+TBsrParameter+" set
    114 co_active="+SqlFormatText("N",GesAct)+""+NL+
    115     "
    116 "+NL+
    117     "where co_node  = "+SqlFormatText(node,GesAct)+"
    118 "+NL+
    119     "  and co_model = "+SqlFormatText(id_model,GesAct)+"
    120 "+NL+
    121     "
    122 "+NL+
    123     "and co_parameter in
    124 "+NL+
    125     "(
    126 "+NL+
    127     "
    128 "+NL+
    129     TxtListItemQuote(lcw,","+NL)+"
    130 "+NL+
    131     "
    132 "+NL+
    133     ")";
    134 
    135     Text WriteFile(logRoute,txt2File);
    136 
    137     Real False
    138   })
    139 };
    140 //////////////////////////////////////////////////////////////////////////////
    141 PutDescription("It takes a set of series and check them to find linear
    142 combinations. If thre are linear combinations, it gives an error message
    143 and
    144 writes a query in a log file in the route logRoute.",
    145 LinComWarning);
    146 //////////////////////////////////////////////////////////////////////////////
    147 
    148 
    149 //////////////////////////////////////////////////////////////////////////////
    150 Set LinComWarningSerie(Set OutputInfo, Set inputs)
    151 //////////////////////////////////////////////////////////////////////////////
    152 {
    153   Set  series   = inputs;
    154   Date TruncIni = OutputInfo::IniEstim;
    155   Date TruncEnd = OutputInfo::EndEstim;
    156 
    157   Polyn dif = SetProd(EvalSet(OutputInfo::Arima,Polyn (Set p){ p[4] }));
    158 
    159   Set seriesTrunc = EvalSet(series,Serie (Serie ser)
    160   {
    161     Serie ser2 = SubSer(ser,TruncIni,TruncEnd);
    162     Serie dif:(ser2/MaxS(Abs(ser2)))
    163   });
    164 
    165   Matrix A = MatSetSeries(seriesTrunc);
    166 
    167   // SVD algorithm builds a set of matrix, the second matrix has the
    168   // eigenvalues.
    169   /* Set svd1 = SVDDecomposition(A*Tra(A)); */ // Luis -> SVD no es
    170 estable
    171   Set svd1 = SVD(A*Tra(A));
    172 
    173   // We transform de second matrix into a rounded matrix because we need
    174   // to search eigenvalues = 0. The 10-10 edge is 1000 times the cholesky
    175   // rounding error so i understand that something below this limit is
    176 negligible
    177   Set svd2 = RoundSVD(svd1,False,0.5,10^(-10));
    178 
    179   // we take de main diagonal of the eigenvalues matrix and we check if
    180 exist
    181   // eigenvalues = 0
    182   Matrix diag = SubDiag( svd2[2] ,0 );
    183   Real   prod = MatSet(diag)[1][Card(MatSet(diag)[1])];
    184 
    185   Set lincom = If(BinEQ(prod,?),
    186   {
    187     Set Lin.Com.Procedure(series,svd2)
    188   },
    189   {
    190     Set If(prod,
    191     {
    192       Set Empty
    193     },
    194     {
    195       Set Lin.Com.Procedure(series,svd2)
    196     })
    197   })
    198 };
    199 //////////////////////////////////////////////////////////////////////////////
    200 PutDescription("It takes a set of series and gives you the number of
    201 linear
    202 combinations that happens",
    203 LinComWarningSerie);
    204 //////////////////////////////////////////////////////////////////////////////
    205 
    206 //////////////////////////////////////////////////////////////////////////////
    207 Set LinComWarningMatrix(Set inputs)
    208 //////////////////////////////////////////////////////////////////////////////
    209 {
    210   Set vectors = inputs;
    211 
    212   Matrix A = Tra(BinGroup("|",vectors));
    213 
    214   // SVD algorithm builds a set of matrix, the second matrix has the
    215   // eigenvalues.
    216   /* Set svd1 = SVDDecomposition(A*Tra(A)); */ // Luis -> SVD no es
    217 estable
    218   Set svd1 = SVD(A*Tra(A));
    219 
    220   // we transform de second matrix into a rounded matrix because we need
    221   // to search eigenvalues = 0
    222   Set svd2 = RoundSVD(svd1,False,0.5,10^(-10));
    223 
    224   // we take de main diagonal of the eigenvalues matrix and we check if
    225 exist
    226   // eigenvalues = 0
    227   Matrix diag = SubDiag( svd2[2] ,0 );
    228   Real   prod = MatSet(diag)[1][Card(MatSet(diag)[1])];
    229 
    230   Set lincom = If(BinEQ(prod,?),
    231   {
    232     Set Lin.Com.Procedure(vectors,svd2)
    233   },
    234   {
    235     Set If(prod,
    236     {
    237       Set Empty
    238     },
    239     {
    240       Set Lin.Com.Procedure(vectors,svd2)
    241     })
    242   })
    243 };
    244 //////////////////////////////////////////////////////////////////////////////
    245 PutDescription("It takes a set of matrix and gives you the number of
    246 linear
    247 combinations that happens",
    248 LinComWarningMatrix);
    249 //////////////////////////////////////////////////////////////////////////////
    250 
    251 //////////////////////////////////////////////////////////////////////////////
    252 Set Lin.Com.Procedure(Set elements, Set svd)
    253 //////////////////////////////////////////////////////////////////////////////
    254 {
    255   // If exist eigenvalues = 0, we rounded to zero the smallest values of
    256 the
    257   // last column of the third SVD matrix, and then, the not zero values
    258 show
    259   // the linear combinations.
    260 
    261   Matrix  lastCol = SubCol(svd[3],[[ Columns(svd[3]) ]]);
    262   Matrix  LogicLC = Not(Not(lastCol));
    263   Real    maxDiag = MaxMatrix(SubDiag(svd[2],0));
    264 
    265   // We take the linear combinatioNs and store them into a set
    266   Set indices = BinGroup("<<",MatSet(LogicLC));
    267 
    268   Real c = 0;
    269 
    270   Set select  = Select(elements,Real (Anything ser)
    271   {
    272     Real c := c+1;
    273     Real If(maxDiag,indices[c],True)
    274   });
    275 
    276   Set EvalSet(select,Text (Anything x)
    277   {
    278     Text Name(x)
    279   })
    280 };
    281 //////////////////////////////////////////////////////////////////////////////
    282 PutDescription("It needs a set of elements (series or matrix) and the
    283 result
    284 of a svd algorithm, and gives you tha distinct variables of elements Set
    285 that
    286 does linear combinations between them",
    287 Lin.Com.Procedure);
    288 //////////////////////////////////////////////////////////////////////////////
    289 
    290 //////////////////////////////////////////////////////////////////////////////
    291 Set RoundSVD(Set svd,Real trasposed, Real minSparse,Real tolerance)
    292 //////////////////////////////////////////////////////////////////////////////
    293 {
    294   Set EvalSet(svd,Matrix (Matrix matrix)
    295   {
    296     VMatrix vmatrix = Mat2VMat(
    297       matrix,     // matrix to transform
    298       trasposed,  // does it trasposed?
    299       minSparse,  // minSparse
    300       tolerance); // tolerance
    301 
    302     Matrix VMat2Mat(vmatrix)
    303   })
    304 };
    305 //////////////////////////////////////////////////////////////////////////////
    306 PutDescription("It rounds an exit of SVD function to avoid negligible
    307 values.
    308 It needs:
    309   matrix     --> matrix to transform
    310   trasposed  --> does it trasposed? (True or False)
    311   minSparse  --> minSparse use 0.5 by default
    312   tolerance  --> round use 10^(-10) recommended
    313 ",
    314 RoundSVD);
    315 //////////////////////////////////////////////////////////////////////////////
    316 
    317 
    318 //////////////////////////////////////////////////////////////////////////////
    319 Set SVDDecomposition(Matrix m)
    320 //////////////////////////////////////////////////////////////////////////////
    321 {
    322   VMatrix vm = Mat2VMat(m);
    323 
    324   Real rows    = VRows(vm);
    325   Real columns = VColumns(vm);
    326   Real nonNull = VNonNullCells(vm);
    327 
    328   Real density = nonNull/(rows*columns);
    329 
    330  // Set If(density <= 0.05, SVD(m,"Sparse"), SVD(m,"Jacobi"))
    331   Set If(density <= 0.01, SVD(m,"Sparse"), SVD(m))
    332 };
    333 //////////////////////////////////////////////////////////////////////////////
    334 PutDescription("It does the SVD decomposition choosing betwwen Sparse and
    335 Jacobi method according to the density of de matrix m.
    336 ",
    337 SVDDecomposition);
    338 //////////////////////////////////////////////////////////////////////////////
    339 
    340 }}}
     6Una cosa a tener en cuenta. El proceso en general funcionaba en torno a la función SVD de tol, sin embargo, como esta función no es exacta sino aproximada, se probaron los diferentes métodos: Golub_Reinsch, Jacobi y Sparse. De todos ellos, el más fiable fué el primero. El Sparse era muy rápido cuando las matrices a factorizar tenían densidades pequeñas (<0.01), y el Jacobi mostraba una salto de nivel muy claro entre los autovalores de las variables que no son combinaciones lineales de los que sí lo son. Sin embargo éstos dos últimos fallaban a veces, detectando combinaciones lineales donde no las había o al revés. Por defecto se
     7decidió usar Golub_Reinsch con un redondeo en 10-10.