00001 #ifndef __TSOPT_EXTOP__
00002 #define __TSOPT_EXTOP__
00003
00004 #include "parserpp.hh"
00005 #include "linop_base.hh"
00006
00007 #ifdef IWAVE_USE_MPI
00008 #include "mpigridpp.hh"
00009 #else
00010 #include "gridpp.hh"
00011 #endif
00012
00013 namespace TSOpt {
00014
00015 using TSOpt::GridSpace;
00016 using namespace RVL;
00017
00018 template<typename Scalar>
00019 class ExtOp: public LinearOp<Scalar> {
00020 private:
00021
00022 GridSpace const & dom;
00023 GridSpace const & rng;
00024 grid const & gdom;
00025 grid const & grng;
00026 string datatype;
00027
00028 FILE * str;
00029
00030
00031 ExtOp();
00032
00033 protected:
00034
00035 void apply(const Vector<Scalar> & x,
00036 Vector<Scalar> & y) const {
00037 try {
00038
00039 PARARRAY * xpar = ps_new();
00040 PARARRAY * ypar = ps_new();
00041 AssignParams xap(*xpar,str);
00042 AssignParams yap(*ypar,str);
00043
00044
00045 x.eval(xap); y.eval(yap);
00046
00047
00048 string xname, yname;
00049 if (!parse<string>(*xpar,datatype,xname)) {
00050 RVLException e;
00051 e<<"Error: ExtOp::apply - failed to extract input filenames\n";
00052 throw e;
00053 }
00054 if (!parse<string>(*ypar,datatype,yname)) {
00055 RVLException e;
00056 e<<"Error: ExtOp::apply - failed to extract output filenames\n";
00057 throw e;
00058 }
00059
00060
00061 string cmd, axisno,smplno;
00062 string tmpf1, tmpf2;
00063 stringstream con1,con2;
00064 tmpf1 = xname; tmpf2 = "spraytmp.rsf";
00065
00066
00067 if(grng.gdim > grng.dim){
00068 con1<<grng.dim+1; axisno = con1.str();
00069 con2<<grng.axes[grng.dim].n; smplno = con2.str();
00070
00071 cmd = "$RSFROOT/bin/sfspray < " + tmpf1 + " > " + tmpf2 + " axis=" + axisno + " n=" + smplno;
00072
00073 if(system(cmd.c_str())){
00074 RVLException e;
00075 e<<"Error: ExtOp::apply - failed to excute sfspray\n";
00076 throw e;
00077 }
00078
00079 tmpf1=tmpf2; tmpf2="spraytmp1.rsf";
00080
00081 for(int i=grng.dim+2;i<=grng.gdim;++i){
00082 con1.str(""); con1<<i; axisno = con1.str();
00083 con2.str(""); con2<<grng.axes[i-1].n; smplno = con2.str();
00084
00085 cmd = "$RSFROOT/bin/sfspray < " + tmpf1 + " > " + tmpf2 + " axis=" + axisno + " n="+ smplno;
00086 tmpf1.swap(tmpf2);
00087 if(system(cmd.c_str())){
00088 RVLException e;
00089 e<<"Error: ExtOp::apply - failed to excute sfspray\n";
00090 throw e;
00091 }
00092 }
00093 con1.str(""); con1<<grng.dim; axisno = con1.str();
00094 con2.str(""); con2<<grng.gdim; smplno = con2.str();
00095
00096
00097 cmd = "$RSFROOT/bin/sfput < " + tmpf1 + " > " + tmpf2 +
00098 " dim=" + axisno + "gdim=" + smplno;
00099 if(system(cmd.c_str())){
00100 RVLException e;
00101 e<<"Error: ExtOp::apply - failed to excute sfput\n";
00102 throw e;
00103 }
00104 }
00105
00106 Vector<Scalar> z(y.getSpace());
00107 AssignFilename af(tmpf2);
00108 z.eval(af);
00109 y.copy(z);
00110
00111
00112 system("/bin/rm spraytmp* $DATAPATH/spraytmp.rsf@");
00113 ps_delete(&xpar); ps_delete(&ypar);
00114 }
00115 catch (RVLException & e) {
00116 e<<"\ncalled in ExtOp::apply\n";
00117 throw e;
00118 }
00119
00120 }
00121
00122
00123 void applyAdj(const Vector<Scalar> & x,
00124 Vector<Scalar> & y) const {
00125 try {
00126
00127
00128 PARARRAY * xpar = ps_new();
00129 AssignParams xap(*xpar,str);
00130
00131
00132 x.eval(xap);
00133
00134
00135 string xname;
00136 string yname;
00137 if (!parse<string>(*xpar,datatype,xname)) {
00138 RVLException e;
00139 e<<"Error: ExtOp::applyAdj - failed to extract filenames\n";
00140 throw e;
00141 }
00142
00143
00144
00145 string cmd;
00146 string axisno, tmpf1, tmpf2;
00147 stringstream convert;
00148 tmpf1 = xname;
00149 tmpf2 = "stacktmp.rsf";
00150
00151 if ( grng.gdim > grng.dim){
00152 convert<<grng.dim+1;
00153 axisno = convert.str();
00154
00155
00156 cmd = "$RSFROOT/bin/sfstack < " + tmpf1 + " > " + tmpf2 +
00157 " axis=" + axisno + " norm=n";
00158
00159 if(system(cmd.c_str())){
00160 RVLException e;
00161 e<<"Error: ExtOp::apply - failed to excute sfstack\n";
00162 throw e;
00163 }
00164 tmpf1=tmpf2;
00165 tmpf2="stacktmp1.rsf";
00166 for(int i=grng.dim+2;i<=grng.gdim;++i){
00167 cmd = "sfstack < " + tmpf1 + " > " + tmpf2 +
00168 " axis=" + axisno + " norm=n";
00169 if(system(cmd.c_str())){
00170 RVLException e;
00171 e<<"Error: ExtOp::apply - failed to excute sfstack\n";
00172 throw e;
00173 }
00174 tmpf1.swap(tmpf2);
00175 }
00176 }
00177
00178 Vector<Scalar> z(y.getSpace());
00179 AssignFilename af(tmpf1);
00180 z.eval(af);
00181 y.copy(z);
00182
00183
00184 system("/bin/rm stacktmp* $DATAPATH/stacktmp.rsf@");
00185 ps_delete(&xpar);
00186 }
00187 catch (RVLException & e) {
00188 e<<"\ncalled in ExtOp::applyAdj\n";
00189 throw e;
00190 }
00191 }
00192
00193 public:
00194
00195 ExtOp(ExtOp const & A):
00196 dom(A.dom), rng(A.rng), str(A.str),
00197 gdom(A.gdom), grng(A.grng),
00198 datatype(A.datatype){
00199 }
00200
00201
00202 ExtOp(GridSpace const & _dom,
00203 GridSpace const & _rng,
00204 string dtype,
00205 FILE * _str):
00206 dom(_dom), rng(_rng), str(_str),
00207 gdom(dom.getGrid()),grng(rng.getGrid()),
00208 datatype(dtype) {
00209
00210 }
00211
00212 ~ExtOp() {}
00213
00214
00215
00216 LinearOp<Scalar> * clone() const { return new ExtOp(*this); }
00217
00218
00219 const Space<float> & getDomain() const { return dom; }
00220 const Space<float> & getRange() const { return rng; }
00221
00222 ostream & write(ostream & str) const {
00223 str<<"test operator for system call op design\n";
00224 return str;
00225 }
00226
00227 };
00228 }
00229 #endif