lsq.c

射线追踪程序

                                        basename=argv[++i];
                                        break;
                                case 'D' :
                                        normdamp=1;
                                case 'd' :
                                        if (argc<=i+1) goto usage;
                                        if (!sscanf(argv[++i],"%f",&damp0)) goto error;
                                        if (damp0<0.) goto error;
                                        break;
                                default  :
                                        goto error;
                        }
                } else goto error;
        }
        if (!basename) goto usage;
        fprintf(stderr,"=================== lsq ===================\n");
        sprintf(sysname,"%s.sys",basename);
        sprintf(paraname,"%s.par",basename);
        sprintf(inmname,"%s.inm",basename);
        sprintf(resolname,"%s.res",basename);
        sprintf(logname,"%s.log",basename);
        sprintf(varname,"%s.var",basename);
        sprintf(outname,"%s.out",basename);

        if (damp0<0.) {
                fprintf(stderr,"Warning: damp=%f\n",damp0=0.);
        }

        if (NULL==(fmat=fopen(sysname,"rb"))) {
                perror(sysname);
                exit(1);
        }
        G=readmatrix(fmat);
        if (!G) {
                exit(1);
        }

        N=G->nl;
        M=G->nc;

        d=(float*)calloc(N,sizeof(float));
        sd=(float*)calloc(N,sizeof(float));
        d2=(float*)calloc(N,sizeof(float));
        m=(float*)calloc(M,sizeof(float));
        dampfact=(float*)calloc(M,sizeof(float));
        res=(float*)calloc(M,sizeof(float));
        sm=(float*)calloc(M,sizeof(float));
        if ((!d)||(!m)||(!sd)||(!sm)) {
                fprintf(stderr,"malloc error (x)\n");
                exit(1);
        }
        fprintf(stderr,"Reading [d]... ");
        if (fread(d,sizeof(float),N,fmat)!=N) {
                fprintf(stderr,"Error (2) reading %s\n",sysname);
                exit(1);
        }
        fprintf(stderr,"[sd]... ");
        if (fread(sd,sizeof(float),N,fmat)!=N) {
                fprintf(stderr,"Error (3) reading %s\n",sysname);
                exit(1);
        }
        fprintf(stderr,"[damp]... ");
        if (fread(dampfact,sizeof(float),M,fmat)!=M) {
                fprintf(stderr,"Error (4) reading %s\n",sysname);
                exit(1);
        }

        fprintf(stderr,"Alloc...");

        Gtd=(float*)calloc(M,sizeof(float));
        GtG=(float**)calloc(M,sizeof(float*));
        inv=(float**)calloc(M,sizeof(float*));
        R=(float**)calloc(M,sizeof(float*));
        for (i=0; i<M; i++) {
                if ((inv[i]=(float*)calloc(i+1,sizeof(float)))==NULL) {
                        fprintf(stderr,"malloc error (2:%d)\n",i);
                        exit(1);
                }
                if ((GtG[i]=(float*)calloc(i+1,sizeof(float)))==NULL) {
                        fprintf(stderr,"malloc error (2:%d)\n",i);
                        exit(1);
                }
        }

        for (i=0; i<M; i++) {
                if ((R[i]=(float*)calloc(M,sizeof(float)))==NULL) {
                        fprintf(stderr,"malloc error (4:%d)\n",i);
                        exit(1);
                }
        }
        fprintf(stderr,"OK\n");

        fprintf(stderr,"Max alloc. mem.: %dMb\n",
                        (
                        sizeof(float)*(
                        4*N + 5*M + M*M + 4*M*(M+1)/2 )+
                        5*M*sizeof(float*)
                        ) / 1024 / 1024);


        fprintf(stderr,"\n\nvar(d)=%10.3e",var=cvar(d,N));

        if (normdamp) for (j=0; j<M; j++) dampfact[j]*=damp0;
        else for (j=0; j<M; j++) dampfact[j]=damp0;


        /*****************************************************/
        fprintf(stderr,"\n   Gt*G, Gt*d... ");
        ptime(1,1);
        ptime(0,1);

        for (k=0; k<N; k++) {
                PROGRESS(k,N);
                for (i=0; i<M; i++) {
                        if (sparsegetval(G,k,i,&x)) {
                                w=GtG[i];
                                for (j=0; j<=i; j++) {
                                        if (sparsegetval(G,k,j,&y)) {
                                                w[j]+=x*y;
                                        }
                                }
                                Gtd[i]+=x*d[k];
                        }
                }
        }

        ptime(0,2);

        /* copie de GtG dans inv car inv_cholesky travaillera "en place" */
        for (i=0; i<M; i++) {
                memcpy(inv[i],GtG[i],(i+1)*sizeof(float));
        }

        fprintf(stderr,"\n   inv(Gt*G+D*Id)... ");
        for (i=0; i<M; i++) {
                inv[i][i]+=dampfact[i];
        }

        fprintf(stderr,"cholesky... ");

        ptime(0,1);
        inv_cholesky(inv,M);
        ptime(0,2);


        ptime(0,1);
        fprintf(stderr,"\n   (inv(Gt*G)+D*Id) * [Gt*d,Gt*G]... ");

        for (k=0; k<M; k++) {
                PROGRESS(k,M);
                v=inv[k];
                for (j=0; j<M; j++) {
                        m[k]+=((k<j)?inv[j][k]:v[j])*Gtd[j];
                        x=((j<k)?GtG[k][j]:GtG[j][k]);
                        if (x>EPS) {
                                for (i=0; i<k; i++) {
                                        R[i][j]+=v[i]*x;
                                }
                                for (; i<M; i++) {
                                        R[i][j]+=inv[i][k]*x;
                                }
                        }
                }
        }

        ptime(0,2);

        ptime(0,1);
        fprintf(stderr,"\n   diag(cov)... ");
        ptime(0,1);
        for (i=0; i<M; i++) sm[i]=0.;
        for (i=0; i<M; i++)
                for (k=0; k<M; k++)
                        sm[i]+=((k<i)?inv[i][k]:inv[k][i])*R[i][k];
        ptime(0,2);

        fprintf(stderr,"\nDone ");
        ptime(1,2);
        for (i=0; i<M; i++) free(inv[i]);
        free(inv);
        for (i=0; i<M; i++) free(GtG[i]);
        free(GtG);

        fprintf(stderr,"\nWriting model (%s)... ",paraname);
        if (NULL==(f=fopen(paraname,"wt"))) {
                perror(paraname);
                exit(1);
        }
        for (j=0; j<M; j++) {
                fprintf(f,"%d %e %e %e\n",j,m[j],sd[i]*sqrt(sm[j]),R[j][j]);
        }
        fclose(f);

        fprintf(stderr," resolution (%s)... ",resolname);

        dump(resolname,R,M,M);

        memset(d2,0,N*sizeof(float));
        for (i=0; i<N; i++) {
                for (j=0; j<M; j++) {
                        if (sparsegetval(G,i,j,&x))
                                d2[i]+=x*m[j];
                }
        }
        fprintf(stderr,"done\n");

        fprintf(stderr,"var(d)=");
        for (i=0; i<N; i++) {
                d2[i]-=d[i];
        }

        fprintf(stderr,"%10.3e\n",var=cvar(d2,N));

        fprintf(stderr,"var(m)=");
        fprintf(stderr,"%10.3e\n",var=cvar(m,M));

        trR=0.;
        for (i=0; i<M; i++) {
                trR+=R[i][i];
        }

        fprintf(stderr,"tr(R)= %10.3e\n",trR);

        fprintf(stderr,"Writing residuals (%s)...\n",outname);
        if (NULL==(f=fopen(outname,"wt"))) {
                perror(outname);
                exit(1);
        }
        for (j=0; j<N; j++) {
                fprintf(f,"%d %e %e\n",j,d[j],d2[j]);
        }
        fclose(f);

        if (VAR) {
                FILE *fvar=fopen(varname,"at");
                float vi=cvar(d,N);
                float vd=cvar(d2,N);
                float vm=cvar(m,M);
                if (!fvar) {
                        perror(varname);
                        exit(0);
                }
                fprintf(fvar,"%f %f %f %f %f %d %d\n",vi,vd,vm,trR,damp0,N,M);
                fclose(fvar);
        }

        return(0);
error:
        fprintf(stderr,"%s: bad argument: '%s'\n",argv[0],argv[i]);
usage:
        fprintf(stderr,"%s -n name [-d|D damp] [-r] [-s] [-v]\n",
                        argv[0]);
        fprintf(stderr,"-D damp=theta2*facteur(name.lay)\n");
        fprintf(stderr,"-d damp=theta2\n");
        fprintf(stderr,"-v MAJ variance (name.var)\n");
        fprintf(stderr,"-s sauvegarde de la matrice inverse\n");
        return(1);
}