diff --git a/meldata.cc b/meldata.cc index 0dc130db0e1665da769bb33fa145bacb425eb82f..02cfe748c80296d5a50c3bf70db7cd89144c2821 100644 --- a/meldata.cc +++ b/meldata.cc @@ -413,16 +413,21 @@ namespace Melodic{ // ICadjust = IC; } else{ + Matrix resids = stdev(Data - mixMatrix * IC); for(int ctr=1;ctr<=resids.Ncols();ctr++) if(resids(1,ctr) < 0.05) resids(1,ctr)=1; - stdNoisei = pow(stdev(Data - mixMatrix * IC)* + // stdNoisei = pow(stdev(Data - mixMatrix * IC)* + // std::sqrt((float)(Data.Nrows()-1))/ + // std::sqrt((float)(Data.Nrows()-IC.Nrows())),-1); + stdNoisei = pow(resids* std::sqrt((float)(Data.Nrows()-1))/ std::sqrt((float)(Data.Nrows()-IC.Nrows())),-1); - + ColumnVector diagvals; diagvals=pow(diag(unmixMatrix*unmixMatrix.t()),-0.5); + save4D(SP(IC,diagvals*stdNoisei),opts.outputfname.value() + "_IC"); } diff --git a/melica.cc b/melica.cc index b5d9b7c8e02105fd7c9553c9c65fd545aafa5281..99568374edcd6be69c82b025bf8d8b38337f7fdd 100644 --- a/melica.cc +++ b/melica.cc @@ -101,9 +101,11 @@ namespace Melodic{ } //termination condition : angle between old and new < epsilon minAbsSin = 1 - diag(abs(redUMM.t()*redUMM_old)).Minimum(); + message(" Step no. " << cum_itt + itt_ctr << " change : " << minAbsSin << endl); - if((abs(minAbsSin) < opts.epsilon.value())&& - (opts.approach.value()!=string("tica"))){ break;} + // if((abs(minAbsSin) < opts.epsilon.value())&& + // (opts.approach.value()!=string("tica"))){ break;} + if((abs(minAbsSin) < opts.epsilon.value())){ break;} itt_ctr++; } while(itt_ctr < opts.maxNumItt.value()); @@ -119,7 +121,7 @@ namespace Melodic{ temp = melodat.reduce_dimred(temp); redUMM = melodat.get_white() * temp; - message(" change : " << abs(minAbsSin2) << endl); + message(" change : " << minAbsSin2 << endl); if(abs(minAbsSin2) < opts.epsilonS.value() && abs(minAbsSin) < opts.epsilon.value()){ break;} } } while( @@ -467,6 +469,7 @@ namespace Melodic{ melodat.set_mix(tmp); melodat.set_IC(temp); + melodat.set_ICstats(scales); melodat.sort(); diff --git a/melodic.cc b/melodic.cc index d6d95d8431bb31bb58af5165b25ad2265e422b0d..0d04f18e8bdbdbd4eb344d13109351de76bef407 100644 --- a/melodic.cc +++ b/melodic.cc @@ -114,7 +114,7 @@ int main(int argc, char *argv[]){ opts.pca_dim.set_T(opts.pca_dim.value()-retry*opts.retrystep); } else{ - if(opts.pca_dim.value()-retry*opts.retrystep < melodat.data_dim()){ + if(opts.pca_dim.value()+retry*opts.retrystep < melodat.data_dim()){ opts.pca_dim.set_T(opts.pca_dim.value()+retry*opts.retrystep); }else{ leaveloop = true; //stupid, but break does not compile @@ -234,7 +234,7 @@ void mmonly(Log& logger, MelodicOptions& opts, message(" done" << endl); } - message("Reading mixing matrix " << opts.filtermix.value()); + message("Reading mixing matrix " << opts.filtermix.value() << " ... "); mixMatrix = read_ascii_matrix(opts.filtermix.value()); if (mixMatrix.Storage()<=0) { cerr <<" Please specify the mixing matrix correctly" << endl; @@ -257,13 +257,12 @@ void mmonly(Log& logger, MelodicOptions& opts, } } - melodat.tempInfo = ICvolInfo; melodat.set_mask(Mask); melodat.set_mean(Mean); melodat.set_IC(ICs); melodat.set_mix(mixMatrix); - fmixMatrix = calc_FFT(melodat.expand_mix(), opts.logPower.value()); + fmixMatrix = calc_FFT(mixMatrix, opts.logPower.value()); melodat.set_fmix(fmixMatrix); fmixMatrix = pinv(mixMatrix); melodat.set_unmix(fmixMatrix); @@ -272,7 +271,7 @@ void mmonly(Log& logger, MelodicOptions& opts, Matrix mmres; Matrix pmaps;//(ICs); - if(opts.perf_mm.value()) + if(opts.perf_mm.value()) mmres = mmall(logger,opts,melodat,report,pmaps); }