Разложение импульса.Разложение ЭДС.Использование регуляризации, страница 15

{

sum = 0;

for( j = 1; j < IM[i+1]-IM[i]; j++ )sum += M[IM[i]+j]*HEMF[i+j];

if( fabs(M[IM[i]]) < maxelem )return errCannotSolveSLAE;

HEMF[i] = (EMF[i] - sum)/M[IM[i]];

}

return errOk;

}

int CEMFTransform::RestHvsLinReg( double *HT, double *HEMF, int &hn )

{

double *betta;

double *A,*F;

double *M;

double mindiff;

double sum, maxelem;

int *IM;

int nx;

int i,j;

int fp, previp, iter;

int res;

res = DelImpulsTime();

if( res != errOk )return res;

if( m < 2 || n < 2 )return errTooFewPoints;

res = 2*m*sizeof(double);

A = (double *)&IntMem[res];

F = (double *)&IntMem[res+=(n+1)*(n+2)/2*sizeof(double)];

betta = (double *)&IntMem[res+=(n+1)*sizeof(double)];

M = (double *)&IntMem[res+=(n+1)*sizeof(double)];

IM = (int *)&IntMem[res+=(n+1)*(n+2)/2*sizeof(double)];

res = MakeMatrix( M, IM, HT, hn );

if( res != errOk )return res;

// Do default solve

maxelem = fabs(M[0]);

for( i = 1; i <= IM[hn-1]; i++ )

maxelem = max( fabs(M[i]), maxelem );

maxelem *= 1e-15;

for( i = hn-1; i >= 0; i-- )

{

sum = 0;

for( j = 1; j < IM[i+1]-IM[i]; j++ )sum += M[IM[i]+j]*HEMF[i+j];

if( fabs(M[IM[i]]) < maxelem )return errCannotSolveSLAE;

HEMF[i] = (EMF[i] - sum)/M[IM[i]];

}

// Check solve

nx = 0; i = hn-1;

while( i >= 0 && HEMF[i] >= HEMF[i+1] ){ nx++; i--; }

if( i > 0 )nx -= (int)(0.05*nx);else return errOk;

mindiff = fabs(HEMF[hn]-HEMF[hn-1])/(HT[hn]-HT[hn-1]);

for( i = hn-nx; i < hn; i++ )

mindiff = min( fabs(HEMF[i]-HEMF[i-1])/(HT[i]-HT[i-1]), mindiff );

// Prepare for regularization

MakeMatrixForReg( A, F, M, IM, hn );

if( ProgressFunc != NULL )

(*ProgressFunc)( paRegular, 0 );

for( i = 0; i < hn-nx; i++ )

betta[i] = bettacoeff*EMF[hn-2]/(HT[i+1]-HT[i]);

for( i = hn-nx; i < hn; i++ )betta[i] = 0;

fp = previp = 0;

iter = 1;

while( fp < hn-nx && iter <= MaxIter )

{

CholeskySR( hn, A, betta, hn-nx, M );

if( !CholeSolveS( hn, M, F, hn-nx, HEMF ) )return errCannotSolveSLAE;

double diff, prevdiff;

prevdiff = fabs(HEMF[1] - HEMF[0])/(HT[1] - HT[0]);

fp = hn-nx;

for( i = 1; i < hn-nx && fp == hn-nx; i++ )

{

diff = fabs(HEMF[i+1] - HEMF[i])/(HT[i+1] - HT[i]);

if( fabs((diff - prevdiff)/max(fabs(prevdiff),mindiff)) > maxdiff )

fp = prevdiff < diff ? i : i-1;else if( (1.-maxinc)*HEMF[i] > HEMF[i-1] )fp = i-1;

prevdiff = diff;

}

if( fp < previp )

for( i = fp; i < previp; i++ )betta[i] *= mulcoeff;else

for( i = fp; i < hn-nx; i++ )betta[i] *= mulcoeff;

previp = fp;

if( ProgressFunc != NULL )

if( !(*ProgressFunc)( paContinue, fp*100/(hn-nx) ) )return errUserBreak;

iter++;

}

if( iter > MaxIter )return errRegFailed;

return errOk;

}

int CEMFTransform::RestHvsLinRegDiff( double *DEMF, double *HT, double *HEMF, int &hn )

{

double *betta;

double *A,*F;

double *M;

double sum, maxelem;

int *IM;

int nx;

int i,j;

int fp, previp, iter;

int res;

res = DelImpulsTime();

if( res != errOk )return res;

if( m < 2 || n < 2 )return errTooFewPoints;

res = 2*m*sizeof(double);

A = (double *)&IntMem[res];

F = (double *)&IntMem[res+=(n+1)*(n+2)/2*sizeof(double)];

betta = (double *)&IntMem[res+=(n+1)*sizeof(double)];

M = (double *)&IntMem[res+=(n+1)*sizeof(double)];

IM = (int *)&IntMem[res+=(n+1)*(n+2)/2*sizeof(double)];

res = MakeMatrix( M, IM, HT, hn );

if( res != errOk )return res;

// Do default solve

maxelem = fabs(M[0]);

for( i = 1; i <= IM[hn-1]; i++ )

maxelem = max( fabs(M[i]), maxelem );

maxelem *= 1e-15;

for( i = hn-1; i >= 0; i-- )

{

sum = 0;

for( j = 1; j < IM[i+1]-IM[i]; j++ )sum += M[IM[i]+j]*HEMF[i+j];

if( fabs(M[IM[i]]) < maxelem )return errCannotSolveSLAE;

HEMF[i] = (EMF[i] - sum)/M[IM[i]];

}

// Check solve

nx = 0; i = hn-1;

while( i >= 0 && HEMF[i] >= HEMF[i+1] ){ nx++; i--; }

if( i > 0 )nx -= (int)(0.05*nx);else return errOk;

// Prepare for regularization

MakeMatrixForReg( A, F, M, IM, hn );

if( ProgressFunc != NULL )

(*ProgressFunc)( paRegular, 0 );