package cquest.preprocessing;

import cquest.quantitation.Svd;
import cquest.Data;
import cquest.rwfile.fileMru;
import java.text.DecimalFormat;
import java.io.*;
import java.util.Vector;
import cquest.utils.MrsiTxtReader;
import cquest.rwfile.FileManagementWrite;

/** class WaterSuppCmdLine.
 * class to suppress the water peak from command line.
 *  
 * @author Michael Sdika, Helene Ratiney 
 */
public class SVDSuppCmdLine
{




 /** Save the new water-suppressed File */

  public void saveResults(Data data, String basename)
    { 
      // save the Signal with watre Suppresion
      double[][][] newSigs = new double[data.getSignalsNb()][][];
      for (int i = 0; i<newSigs.length ; i++) 
	{ newSigs[i] = data.getSignal(i);
       	}

      try
	{ FileManagementWrite.SaveMruFromRawSignalS(
	      basename + "_SVDSup",
	      data.stepTime,
	      data.B0,
	      data.transmitterFreq,
	      data.nucleus,
	      data.hzref,
	      data.ppmref,
	      data.zeroOrderPhase,
	      data.beginTime,
	      data.getOriginalLength(),
	      newSigs);
	}
      catch (IOException e)
	{ System.out.println ("unable to save the new file");
	}
    }    

  
  /** print the usage on stdout. */
  public void printUsage()
    { System.out.println ("usage : waterSupp [-b bounds] [-m mrsi.txt] data1.mrui ...");
      System.out.println ("-m mrsi.txt  : if specified, give the selection of voxels to quantitate");
      System.out.println ("-k nbComponent : number of Component used in HLSVD" );
      System.out.println ("-b bounds1 bounds2 : frequency bounds of the water (in Hz) with bounds1<bounds2" );
      System.out.println ("the data is a .mrui file");
      System.out.println ("option must be given before data files.");
    }

  public void launch(String datafile, String mrsiFile, int k, double[] bounds)    
    {
      Data data = null;

      String outfile = null;
            
      try
      { System.out.println ("datafile="+datafile);
	  outfile = datafile.substring(0, datafile.length()-5);
       System.out.println ("outfile="+outfile);
      }
      catch (IndexOutOfBoundsException e)                            
      { System.out.println (e.getMessage());
        e.printStackTrace();
      return ;
      }

      try
	{ data = (new fileMru(datafile)).loadMru();
          DecimalFormat extensionNumberFormat = new DecimalFormat("0000");
	  for (int i=0;i<data.getSignalsNb();i++)
	    { 
	      data.getFIDAt(i).SetName((new File (outfile)).getName()+"_"+ extensionNumberFormat.format(i)+".mrui");
	    }
	  
	}
      catch (FileNotFoundException e)
	{ System.out.println ("File "+datafile+" not found");
	}  
      catch(Exception e)
	{ System.out.println ("Error when loading file "+datafile);
	  System.out.println (e.getMessage());
	  return ;
	}
        
        boolean[][] mrsiSel = null;
        if (mrsiFile!=null)
	
	{ try
	    { mrsiSel = MrsiTxtReader.read(new File(mrsiFile));
	      data.dataWidthForMrsi  = mrsiSel.length;
	      data.dataHeightForMrsi = mrsiSel[0].length;
	      
	    }
	  catch (Exception e)
	    { System.out.println ("Error while reading mrsi file");
	      System.out.println (e.getMessage());
	      System.exit(1);
	    }
	}

        //new Svd(Svd.HLSVD,data,data.getOriginalLength()/2,k,data.getOriginalLength(),true,false,null, 0);
        Data dataWS=processSVD(data,k,bounds , mrsiSel);
        saveResults(dataWS,outfile);

    }



    protected Data processSVD(Data data, int k, double[] bounds, boolean [][] mrsiSelectionArray) {
 
    Data dataWaterSupp=new Data();
       	
	int iEchoMax = -1;
    if (data.areEchoes) {
      iEchoMax = data.getMaxIndexFromSerie();
    }
    if ((iEchoMax < 2) || (iEchoMax > data.getOriginalLength() - 3)) { // -> FIDs (-1) or invalid echoes
      
        if (data.areEchoes) System.out.println ("Warning: Invalid Echoes... Quantifying as FIDs");
        
        for (int iSignal = 0; iSignal < data.getSignalsNb(); iSignal++) {
	 System.out.println ("Signal "+(iSignal+1) +" out of "+data.getSignalsNb()+" is being processed");	       
          double[][] signalTemp = data.getSignal(iSignal);
          Data dataTemp = new Data();
	  dataWaterSupp.setStepTime(data.stepTime);
          dataWaterSupp.B0=data.B0;
	  dataWaterSupp.setTransmitterFrequency(data.getTransmitterFrequency());
	  dataWaterSupp.nucleus =data.nucleus;
	  dataWaterSupp.setErfrequency(data.getErfrequency());
	  dataTemp.setStepTime(data.stepTime);
	  dataWaterSupp.addSignal(signalTemp);
          dataTemp.addSignal(signalTemp);
	  dataTemp.findResult();
         int dataWidth=1;
	 int dataHeight=1;
	  if (mrsiSelectionArray!=null) {
          dataWidth = data.dataWidthForMrsi;
	  dataHeight= data.dataHeightForMrsi;
	  } 
          
	 if (mrsiSelectionArray!=null&&!mrsiSelectionArray[iSignal % dataWidth][iSignal / dataWidth]) 
	   {	      
            dataTemp.result.addFID(0, data.stepTime);
            dataTemp.result.setResultAt(0, new double[k], new double[k], (double[]) null, new double[k], new double[k]);
            dataTemp.result.setSingularValuesAt(0, new double[k]);	    
           }
	 else 
	   {
	    
	    new Svd(Svd.HLSVD,dataTemp,data.getOriginalLength()/2,k,data.getOriginalLength(),true,false,null, 0);
            int[] numToRemove = new int[dataTemp.result.getNBPeakOfFIDAt(0)];
            int nbOfComp = 0;
            int indice = 0;
            double[] amplSVD = new double[dataTemp.result.getNBPeakOfFIDAt(0)];
            double[] dampSVD = new double[dataTemp.result.getNBPeakOfFIDAt(0)];
            double[] dampRSVD = new double[dataTemp.result.getNBPeakOfFIDAt(0)];
            double[] freqSVD = new double[dataTemp.result.getNBPeakOfFIDAt(0)];
            double[] phasSVD = new double[dataTemp.result.getNBPeakOfFIDAt(0)];
            double[] dampSVDcramer = new double[dataTemp.result.getNBPeakOfFIDAt(0)];
            double[] freqSVDcramer = new double[dataTemp.result.getNBPeakOfFIDAt(0)];
            double[] phasSVDcramer = new double[dataTemp.result.getNBPeakOfFIDAt(0)];
	    Vector tempF = dataTemp.result.getFrequencyAt(0);
            Vector tempA = dataTemp.result.getAmplitudeAt(0);
            Vector tempD = dataTemp.result.getDampingAt(0);
      	    Vector tempP = dataTemp.result.getPhaseAt(0);
      	    Vector tempAe = dataTemp.result.getAmplitudeErrorAt(0);
      	    Vector tempDe = dataTemp.result.getDampingErrorAt(0);
      	    Vector tempFe = dataTemp.result.getFrequencyErrorAt(0);
      	    Vector tempPe = dataTemp.result.getPhaseErrorAt(0);
            for (int iComponent = 0; iComponent < amplSVD.length; iComponent++) {
              amplSVD[iComponent] = ((Double) tempA.elementAt(iComponent)).doubleValue();
              dampSVD[iComponent] = ((Double) tempD.elementAt(iComponent)).doubleValue() * 1000 / data.stepTime;
      	      dampRSVD[iComponent] = 0.0;
      	      freqSVD[iComponent] = ((Double) tempF.elementAt(iComponent)).doubleValue() * 1000 / data.stepTime;
      	      phasSVD[iComponent] = ((Double) tempP.elementAt(iComponent)).doubleValue() * 180 / Math.PI;       
	      numToRemove[iComponent] = 0;
	     for (int l=0;l<bounds.length;l=l+2)
	       { 
	       if (freqSVD[iComponent] >= bounds[l] && freqSVD[iComponent] <= bounds[l+1])		 
                 {
                 numToRemove[iComponent] = 1;
                 nbOfComp ++;
                 }	    
	       }
	   }
           
	   if (nbOfComp > 0) 
	     {	       
	      for (int iComponent = 0; iComponent < numToRemove.length; iComponent++) 
		{
                 if (numToRemove[iComponent] == 1) 
		   {
		    double[][] daComponent = dataTemp.result.getSignalComponent(0)[iComponent];
		     for (int iDataPoint = 0; iDataPoint < dataTemp.getOriginalLength(); iDataPoint++) {
	      		signalTemp[0][iDataPoint] -= daComponent[0][iDataPoint];
	      		signalTemp[1][iDataPoint] -= daComponent[1][iDataPoint];
 		     }
   		   }
  		}
              dataWaterSupp.setDataOfSignal(iSignal, signalTemp);

	     }
	   
         }//end if selected
      }//end for all signal 
    }
    else {//**************** The echo case remain to be treated**************************************//
    System.out.println ("The signal is considered as an echo; for th moment nothing is done in this case !!!");	
    }
    return (dataWaterSupp);
  }


  

  /** main function : to run the HLSVD water suppression from the command line. */

  public static void main(String[] args)
    { 

      SVDSuppCmdLine svdscl = new SVDSuppCmdLine();
      
      if (args.length==0)
	   { svdscl.printUsage();
	      System.exit(1); 
	    }
      String mrsiFile  = null;
      
      int k=25;
      Vector vectBounds = new Vector();
      //vectBounds.addElement( new Double(-20));
      //vectBounds.addElement( new Double(20));
      
      
      int next = 0;
      boolean optiongiven = false;
      while (next<args.length && !optiongiven)
      { if (args[next].equals("-m"))
	  { if (next+1==args.length)
		{ System.out.println ("mrsi file not given");
	          svdscl.printUsage();
		  System.exit(1);
		  
		}
	      mrsiFile  = args[next+1];
	      next += 2;
	  }
  	  else if (args[next].equals("-k"))
	    { if (next+1==args.length)
		{ System.out.println ("number of Components not given");
	          svdscl.printUsage();
		  System.exit(1);
		  
		}
	      k  = (new Integer(args[next+1])).intValue();
	      next += 2;
	    }
	  
	  else if (args[next].equals("-b"))
	    { if (next+2>=args.length)
		{ System.out.println ("frequency bounds  not given");
	          svdscl.printUsage();
		  System.exit(1);
		  
		}
              double arg0= new Double((args[next+1])).doubleValue(); 
	      double arg1= new Double((args[next+2])).doubleValue();     
	      vectBounds.addElement(new Double(arg0));
              vectBounds.addElement(new Double(arg1));

	      if (arg0>=arg1)
		{ System.out.println ("frequency bounds  not given in the good order b1<b2 or are equal");
	          svdscl.printUsage();
		  System.exit(1);
		}
	      

	        
	      	      
	      next += 3;
	    }
	  else
	    { optiongiven = true;
	    }
      }	


      double [] bounds=new double [vectBounds.size()]; 
      for (int l=0;l<vectBounds.size();l=l+2)
	{ 
	bounds[l] = ((Double)(vectBounds.elementAt(l))).doubleValue();
	bounds[l+1] = ((Double)(vectBounds.elementAt(l+1))).doubleValue();
      	
      	System.out.println ("bounds "+bounds[l]+" "+bounds[l+1]);
	}
      
      for (int i=next;i<args.length;i++)
	{svdscl.launch(args[i],mrsiFile,k,bounds);
	}
  
  
  /*Data d, int r, int k, int points, boolean all, boolean calculateErrors, boolean isLanczos, boolean isTotalLS, int noiseFromResidu) 
    int method;
    if (isLanczos)
      { method = isTotalLS ? Svd.HLTLS : Svd.HLSVD ;
      }
    else
      { method = isTotalLS ? Svd.HTLS  : Svd.HSVD ;
      }
    
    new Svd(method,data,r,k,points,all,calculateErrors,null, 0);*/
  }


}