/* var_test.c: 
 
   Compilation from Matlab:
   mex object_update_norm.c
   maybe a tiny bit faster code is generated by 
   mex -O COPTIMFLAGS='-O2' LDOPTIMFLAGS='-O2' object_update_norm.c
 
   Usage from Matlab:
   object_update_norm(iter,cprobe,ob,pr_nrm,int32(positions),int32(numpts));
 
 This code in matlab:
 asize = size(cprobe);
 for i=1:numpts
    Indy = positions(i,1) + (1:asize(1));
    Indx = positions(i,2) + (1:asize(2));
    ob(Indy,Indx) = ob(Indy,Indx) + cprobe .* iter(:,:,i);
    pr_nrm(Indy,Indx) = pr_nrm(Indy,Indx) + abs(cprobe).^2;
 end
 
Academic License Agreement

Source Code

 Introduction 
 •	This license agreement sets forth the terms and conditions under which the PAUL SCHERRER INSTITUT (PSI), CH-5232 Villigen-PSI, Switzerland (hereafter "LICENSOR") 
   will grant you (hereafter "LICENSEE") a royalty-free, non-exclusive license for academic, non-commercial purposes only (hereafter "LICENSE") to use the cSAXS 
   ptychography MATLAB package computer software program and associated documentation furnished hereunder (hereafter "PROGRAM").

 Terms and Conditions of the LICENSE
 1.	LICENSOR grants to LICENSEE a royalty-free, non-exclusive license to use the PROGRAM for academic, non-commercial purposes, upon the terms and conditions 
       hereinafter set out and until termination of this license as set forth below.
 2.	LICENSEE acknowledges that the PROGRAM is a research tool still in the development stage. The PROGRAM is provided without any related services, improvements 
       or warranties from LICENSOR and that the LICENSE is entered into in order to enable others to utilize the PROGRAM in their academic activities. It is the 
       LICENSEE’s responsibility to ensure its proper use and the correctness of the results.”
 3.	THE PROGRAM IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR 
       A PARTICULAR PURPOSE AND NONINFRINGEMENT OF ANY PATENTS, COPYRIGHTS, TRADEMARKS OR OTHER RIGHTS. IN NO EVENT SHALL THE LICENSOR, THE AUTHORS OR THE COPYRIGHT 
       HOLDERS BE LIABLE FOR ANY CLAIM, DIRECT, INDIRECT OR CONSEQUENTIAL DAMAGES OR OTHER LIABILITY ARISING FROM, OUT OF OR IN CONNECTION WITH THE PROGRAM OR THE USE 
       OF THE PROGRAM OR OTHER DEALINGS IN THE PROGRAM.
 4.	LICENSEE agrees that it will use the PROGRAM and any modifications, improvements, or derivatives of PROGRAM that LICENSEE may create (collectively, 
       "IMPROVEMENTS") solely for academic, non-commercial purposes and that any copy of PROGRAM or derivatives thereof shall be distributed only under the same 
       license as PROGRAM. The terms "academic, non-commercial", as used in this Agreement, mean academic or other scholarly research which (a) is not undertaken for 
       profit, or (b) is not intended to produce works, services, or data for commercial use, or (c) is neither conducted, nor funded, by a person or an entity engaged 
       in the commercial use, application or exploitation of works similar to the PROGRAM.
 5.	LICENSEE agrees that it shall make the following acknowledgement in any publication resulting from the use of the PROGRAM or any translation of the code into 
       another computing language:
       "Data processing was carried out using the cSAXS ptychography MATLAB package developed by the Science IT and the coherent X-ray scattering (CXS) groups, Paul 
       Scherrer Institut, Switzerland."

 Additionally, any publication using the package, or any translation of the code into another computing language should cite for difference map:
 P. Thibault, M. Dierolf, A. Menzel, O. Bunk, C. David, F. Pfeiffer, High-resolution scanning X-ray diffraction microscopy, Science 321, 379–382 (2008). 
   (doi: 10.1126/science.1158573),
 for maximum likelihood:
 P. Thibault and M. Guizar-Sicairos, Maximum-likelihood refinement for coherent diffractive imaging, New J. Phys. 14, 063004 (2012). 
   (doi: 10.1088/1367-2630/14/6/063004),
 for mixed coherent modes:
 P. Thibault and A. Menzel, Reconstructing state mixtures from diffraction measurements, Nature 494, 68–71 (2013). (doi: 10.1038/nature11806),
 and/or for multislice:
 E. H. R. Tsai, I. Usov, A. Diaz, A. Menzel, and M. Guizar-Sicairos, X-ray ptychography with extended depth of field, Opt. Express 24, 29089–29108 (2016). 
   (doi: 10.1364/OE.24.029089).
 6.	Except for the above-mentioned acknowledgment, LICENSEE shall not use the PROGRAM title or the names or logos of LICENSOR, nor any adaptation thereof, nor the 
       names of any of its employees or laboratories, in any advertising, promotional or sales material without prior written consent obtained from LICENSOR in each case.
 7.	Ownership of all rights, including copyright in the PROGRAM and in any material associated therewith, shall at all times remain with LICENSOR, and LICENSEE 
       agrees to preserve same. LICENSEE agrees not to use any portion of the PROGRAM or of any IMPROVEMENTS in any machine-readable form outside the PROGRAM, nor to 
       make any copies except for its internal use, without prior written consent of LICENSOR. LICENSEE agrees to place the following copyright notice on any such copies: 
       © All rights reserved. PAUL SCHERRER INSTITUT, Switzerland, Laboratory for Macromolecules and Bioimaging, 2017. 
 8.	The LICENSE shall not be construed to confer any rights upon LICENSEE by implication or otherwise except as specifically set forth herein.
 9.	DISCLAIMER: LICENSEE shall be aware that Phase Focus Limited of Sheffield, UK has an international portfolio of patents and pending applications which relate 
       to ptychography and that the PROGRAM may be capable of being used in circumstances which may fall within the claims of one or more of the Phase Focus patents, 
       in particular of patent with international application number PCT/GB2005/001464. The LICENSOR explicitly declares not to indemnify the users of the software 
       in case Phase Focus or any other third party will open a legal action against the LICENSEE due to the use of the program.
 10.	This Agreement shall be governed by the material laws of Switzerland and any dispute arising out of this Agreement or use of the PROGRAM shall be brought before 
       the courts of Zürich, Switzerland. 
 */

#include "mex.h"
#include <math.h> 
#include <stdio.h> 


void mexFunction(int nlhs, mxArray *plhs[],
                 int nrhs, const mxArray *prhs[])
{
  int i;
    
    /* Check for proper number of arguments. */
  if (nrhs != 6) 
    mexErrMsgIdAndTxt("MexError:ptycho","Six input arguments required: object_update(iter,cprobe,ob,positions,asize,numpts)");
  else if (nlhs != 0) 
    mexErrMsgIdAndTxt("MexError:ptycho","No output argument has to be specified.");
  
    /* Input must be of type double. */
  for (i=0; i < 4; i++) {
      if (mxIsDouble(prhs[i]) != 1){
          printf("Input %d is not double\n",i+1);
          mexErrMsgIdAndTxt("MexError:ptycho","Inputs must be of correct type.");
      }
  }
    /* Input must be of type int32. */
  for (i=4; i<nrhs; i++){
    if (mxIsInt32(prhs[i]) != 1){
          printf("Input %d is not integer\n",i+1);
          mexErrMsgIdAndTxt("MexError:ptycho","Inputs must be of correct type.");
    }
  }
 
    /* It cannot be one-dimensional */
  if(mxGetNumberOfDimensions(prhs[0]) < 2) {
    printf("The 1st input argument must have at least two dimensions.");
    mexErrMsgIdAndTxt("MexError:ptycho","wrong number of dimensions");
  }
    /* It cannot be more than 3-dimensional */
   if(mxGetNumberOfDimensions(prhs[0]) > 3) {
    printf("The 1st input argument must have at most three dimensions.");
    mexErrMsgIdAndTxt("MexError:ptycho","wrong number of dimensions");
  }
  
    /* allocate memory for image data, to be returned 
  plhs[0] = 
    mxCreateNumericMatrix(dim1, dim2, mxDOUBLE_CLASS, mxREAL);
  if (plhs[0] == NULL)
    mexErrMsgIdAndTxt("MexError:ptycho","Could not allocate memory for return data.");*/
  
  double *iter_real, *iter_imag, *cprobe_real, *cprobe_imag, *ob_real, *ob_imag, *pr_nrm;
  int *positions, *numpts;
      
     /* Check that arrays are complex */
  if(mxIsComplex(prhs[0]) != 1) {
      printf("iter input argument must be complex-valued.");
      mexErrMsgIdAndTxt("MexError:ptycho","Expected complex arrays");
  }
  if(mxIsComplex(prhs[1]) != 1) {
      printf("cprobe input argument must be complex-valued.");
      mexErrMsgIdAndTxt("MexError:ptycho","Expected complex arrays");
  }  
  if(mxIsComplex(prhs[2]) != 1) {
      printf("object input argument must be complex-valued.");
      mexErrMsgIdAndTxt("MexError:ptycho","Expected complex arrays");
  }  
     /* get pointers to input data */
  iter_real = mxGetPr(prhs[0]);
  iter_imag = mxGetPi(prhs[0]);
  cprobe_real = mxGetPr(prhs[1]);
  cprobe_imag = mxGetPi(prhs[1]);  
  ob_real = mxGetPr(prhs[2]);
  ob_imag = mxGetPi(prhs[2]); 
  pr_nrm = mxGetPr(prhs[3]);
  positions = mxGetData(prhs[4]);
  numpts = mxGetData(prhs[5]);
  
  /* Get dimension of probe and object */
  int no_rows_p, no_cols_p, no_rows_o, no_cols_o;
  no_rows_p = mxGetM(prhs[1]);
  no_cols_p = mxGetN(prhs[1]);
  no_rows_o = mxGetM(prhs[2]);
  no_cols_o = mxGetN(prhs[2]);
  
  /* Check that iter agrees with probe dimensions and numpts */
  if(mxGetM(prhs[0]) != no_rows_p) {
      mexErrMsgIdAndTxt("MexError:ptycho","iter and cprobe do not have the same number of rows");
  } 
  if(mxGetN(prhs[0]) != numpts[0]*mxGetN(prhs[1])) {
      mexErrMsgIdAndTxt("MexError:ptycho","iter size does not agree with columns in cprobe and numpts");
  }   
  
    /* Check that pr_nrm agrees with ob dimensions */
  if(mxGetM(prhs[3]) != no_rows_o) {
      mexErrMsgIdAndTxt("MexError:ptycho","pr_nrm and ob do not have the same number of rows");
  }
  if(mxGetN(prhs[3]) != no_cols_o) {
      mexErrMsgIdAndTxt("MexError:ptycho","pr_nrm and ob do not have the same number of cols");
  }
  
  /* Check that positions has the right size */
  if(mxGetN(prhs[4]) != 2) {
      mexErrMsgIdAndTxt("MexError:ptycho","positions must have 2 columns");
  } 
  if(mxGetM(prhs[4]) != numpts[0]) {
      mexErrMsgIdAndTxt("MexError:ptycho","positions must have numpts of rows");
  }  
  
  /* Find maximum and minimum of positions */
  int min_pos_row = positions[0], max_pos_row = positions[0];
  int min_pos_col = positions[numpts[0]], max_pos_col = positions[numpts[0]];
  for (i = 1; i < numpts[0] ; i++) {
      if (min_pos_row > positions[i]) {
          min_pos_row = positions[i];
      }
      if (max_pos_row < positions[i]) {
          max_pos_row = positions[i];
      }     
      if (min_pos_col > positions[i+numpts[0]]) {
          min_pos_col = positions[i+numpts[0]];
      }
      if (max_pos_col < positions[i+numpts[0]]) {
          max_pos_col = positions[i+numpts[0]];
      }         
  }
  /*printf("%d min_pos_row\n",min_pos_row);
  printf("%d max_pos_row\n",max_pos_row);
  printf("%d min_pos_col\n",min_pos_col);
  printf("%d max_pos_col\n",max_pos_col);  */
  
  /* Check that indices will not exceed matrix dimensions */
  if ((min_pos_col < 0)||(min_pos_row < 0))  {
      mexErrMsgIdAndTxt("MexError:ptycho","Subscript indices must either be real positive integers\n");
  }
  if ((max_pos_col + no_cols_p > no_cols_o)||(max_pos_row + no_rows_p > no_rows_o))  {
      mexErrMsgIdAndTxt("MexError:ptycho","Subscript exceeds matrix dimensions\n");
  }  
  
  /* The actual code */
  int col=0, row=0, view=0;
  double a2probe [no_rows_p*no_cols_p];
  
  for (col=0;col<no_cols_p;col++) {
      for (row=0;row<no_rows_p;row++) {
          a2probe[row + col*no_rows_p] = pow(cprobe_real[row + col*no_rows_p],2) + pow(cprobe_imag[row + col*no_rows_p],2);
      }
  }
  
  for (view=0;view<numpts[0];view++) {
      for (col=0;col<no_cols_p;col++) {
          for (row=0;row<no_rows_p;row++) {
              /*ob_real[row + positions[view*2] + (col+positions[view*2+1])*no_rows_o] += cprobe_real[row + col*no_rows_p]*iter_real[row + col*no_rows_p + view*no_cols_p*no_rows_p];*/
              ob_real[row + positions[view] + (col+positions[view+numpts[0]])*no_rows_o] += cprobe_real[row + col*no_rows_p]*iter_real[row + col*no_rows_p + view*no_cols_p*no_rows_p] - cprobe_imag[row + col*no_rows_p]*iter_imag[row + col*no_rows_p + view*no_cols_p*no_rows_p];
              ob_imag[row + positions[view] + (col+positions[view+numpts[0]])*no_rows_o] += cprobe_real[row + col*no_rows_p]*iter_imag[row + col*no_rows_p + view*no_cols_p*no_rows_p] + cprobe_imag[row + col*no_rows_p]*iter_real[row + col*no_rows_p + view*no_cols_p*no_rows_p];
              pr_nrm[row + positions[view] + (col+positions[view+numpts[0]])*no_rows_o] += a2probe[row + col*no_rows_p];
              /*printf("%d\n",positions[view*2+1]);*/
          }
      }
  }
  
  return;
}
