rho = scalar. rho_i = 1/(y_i \dot s_i)
184085
@Override
    public void postStep(INDArray gradient) {
        INDArray previousParameters = (INDArray) searchState.get("oldparams");
        INDArray parameters = model.params();
        INDArray previousGradient = (INDArray) searchState.get(GRADIENT_KEY);

        LinkedList<Double> rho = (LinkedList<Double>) searchState.get("rho");
        LinkedList<INDArray> s = (LinkedList<INDArray>) searchState.get("s");
        LinkedList<INDArray> y = (LinkedList<INDArray>) searchState.get("y");

        double sy = Nd4j.getBlasWrapper().dot(previousParameters, previousGradient) + Nd4j.EPS_THRESHOLD;
        double yy = Nd4j.getBlasWrapper().dot(previousGradient, previousGradient) + Nd4j.EPS_THRESHOLD;

        INDArray sCurrent;
        INDArray yCurrent;
        if (s.size() >= m) {
            //Optimization: Remove old (no longer needed) INDArrays, and use assign for re-use.
            //Better to do this: fewer objects created -> less memory overall + less garbage collection
            sCurrent = s.removeLast();
            yCurrent = y.removeLast();
            rho.removeLast();
            sCurrent.assign(parameters).subi(previousParameters);
            yCurrent.assign(gradient).subi(previousGradient);
        } else {
            //First few iterations. Need to allocate new INDArrays for storage (via copy operation sub)
            sCurrent = parameters.sub(previousParameters);
            yCurrent = gradient.sub(previousGradient);
        }

        rho.addFirst(1.0 / sy); //Most recent first
        s.addFirst(sCurrent); //Most recent first. si = currParams - oldParams
        y.addFirst(yCurrent); //Most recent first. yi = currGradient - oldGradient

        //assert (s.size()==y.size()) : "Gradient and parameter sizes are not equal";
        if (s.size() != y.size())
            throw new IllegalStateException("Gradient and parameter sizes are not equal");

        //In general: have m elements in s,y,rho.
        //But for first few iterations, have less.
        int numVectors = Math.min(m, s.size());

        double[] alpha = new double[numVectors];

        // First work backwards, from the most recent difference vectors
        Iterator<INDArray> sIter = s.iterator();
        Iterator<INDArray> yIter = y.iterator();
        Iterator<Double> rhoIter = rho.iterator();

        //searchDir: first used as equivalent to q as per N&W, then later used as r as per N&W.
        //Re-using existing array for performance reasons
        INDArray searchDir = (INDArray) searchState.get(SEARCH_DIR);
        searchDir.assign(gradient);

        for (int i = 0; i < numVectors; i++) {
            INDArray si = sIter.next();
            INDArray yi = yIter.next();
            double rhoi = rhoIter.next();

            if (si.length() != searchDir.length())
                throw new IllegalStateException("Gradients and parameters length not equal");

            alpha[i] = rhoi * Nd4j.getBlasWrapper().dot(si, searchDir);
            Nd4j.getBlasWrapper().level1().axpy(searchDir.length(), -alpha[i], yi, searchDir); //q = q-alpha[i]*yi
        }

        //Use Hessian approximation initialization scheme
        //searchDir = H0*q = (gamma*IdentityMatrix)*q = gamma*q
        double gamma = sy / yy;
        searchDir.muli(gamma);

        //Reverse iterators: end to start. Java LinkedLists are doubly-linked,
        // so still O(1) for reverse iteration operations.
        sIter = s.descendingIterator();
        yIter = y.descendingIterator();
        rhoIter = rho.descendingIterator();
        for (int i = 0; i < numVectors; i++) {
            INDArray si = sIter.next();
            INDArray yi = yIter.next();
            double rhoi = rhoIter.next();

            double beta = rhoi * Nd4j.getBlasWrapper().dot(yi, searchDir); //beta = rho_i * y_i^T * r
            //r = r + s_i * (alpha_i - beta)
            Nd4j.getBlasWrapper().level1().axpy(gradient.length(), alpha[i] - beta, si, searchDir);
        }

        previousParameters.assign(parameters);
        previousGradient.assign(gradient); //Update gradient. Still in searchState map keyed by GRADIENT_KEY
    }
743362
def gaussian2d(xy, sx, sy, mx=0, my=0, rho=0, amp=1, offs=0):
    '''
    
    '''
    x,y = xy
    return offs+amp*( 
        1/(2*np.pi*sx*sy*(1-(rho**2))**0.5) *
         np.exp( (-1/(2*(1-rho**2))) *
                 (
                    ( (x-mx)**2/sx**2 )
                  + ( (y-my)**2/sy**2 )
                  - ( ( 2*rho*(x-mx)*(y-my)) / (sx*sy) )
                  )
                )
         )
1089715
def ystep(self):
        
        """

        rho = self.Nb * self.rho
        mAXU = np.mean(self.AX + self.U, axis=-1)
        self.Y[:] = self.prox_g(mAXU, rho)
959573
def mass_3d(self, r, rho0, Ra, Rs):
        
        m_3d = 4 * np.pi * rho0 * Ra ** 2 * Rs ** 2 / (Rs ** 2 - Ra ** 2) * (Rs * np.arctan(r / Rs) - Ra * np.arctan(r / Ra))
        return m_3d
171276
def _h(self, X, tau):

        

        def cos_func(y):
            if isinstance(y, float) or isinstance(y, int):
                if y > 1:
                    return np.arccosh(y)
                else:
                    return np.arccos(y)
            else:
                values = np.ones_like(y)
                inds1 = np.where(y < 1)
                inds2 = np.where(y > 1)
                values[inds1] = np.arccos(y[inds1])
                # values[inds2] = np.arccosh(y[inds2])
                values[inds2] = np.arccos((1 - y[inds2]) ** .5)
                return values

        def nfw_func(x):
            if isinstance(x,float) or isinstance(x,int):
                if x<1:
                    return np.log(x / 2.) ** 2 - np.arccosh(1. / x) ** 2
                else:
                    return np.log(x / 2.)**2 - np.arccos(1. /x)**2
            else:
                inds1 = np.where(x < 1)
                inds2 = np.where(x >= 1)
                y = np.zeros_like(x)
                y[inds1] = np.log(x[inds1] / 2.) ** 2 - np.arccosh(1. / x[inds1]) ** 2
                y[inds2] = np.log(x[inds2] / 2.) ** 2 - np.arccosh(1. / x[inds2]) ** 2

                return y

        def tnfw_func(x,tau):

            u = x**2
            t2 = tau**2
            Lx = self.L(u**.5, tau)

            return (t2 + 1) ** -2 * (
                2 * t2 * np.pi * (tau - (t2 + u) ** .5 + tau * np.log(tau + (t2 + u) ** .5))
                +
                2 * (t2 - 1) * tau * (t2 + u) ** .5 * Lx
                +
                t2 * (t2 - 1) * Lx ** 2
                +
                4 * t2 * (u - 1) * self.F(u**.5)
                +
                t2 * (t2 - 1) * (cos_func(u ** -0.5)) ** 2
                +
                t2 * ((t2 - 1) * np.log(tau) - t2 - 1) * np.log(u)
                -
                t2 * (
                (t2 - 1) * np.log(tau) * np.log(4 * tau) + 2 * np.log(0.5 * tau) - 2 * tau * (tau - np.pi) * np.log(
                    tau * 2)))

        c = 1e-9
        rescale = 1

        if np.any(X-c<1):

            warnings.warn('Truncated NFW potential not yet implemented for x<1. Using the expression for the NFW '
                          'potential in this regime isntead.')
            rescale = tnfw_func(1.00001,tau)*nfw_func(1)**-1

        if isinstance(X,float) or isinstance(X,int):
            if X<1:
                X = min(X,0.001)
                return nfw_func(X)*rescale
            else:
                return tnfw_func(X,tau)
        else:

            X[np.where(X<0.001)] = 0.001

            values = np.zeros_like(X)
            inds1,inds2 = np.where(X<1),np.where(X>=1)
            values[inds1] = nfw_func(X[inds1])*rescale
            values[inds2] = tnfw_func(X[inds2],tau)

            return values
1327673
def ystep(self):
        
        """

        self.Y = np.asarray(sp.prox_l2(
            self.AX + self.U, (self.lmbda/self.rho)*self.Wtvna,
            axis=self.saxes), dtype=self.dtype)
1084856
def _evaluate(self,R,phi=0.,t=0.):
        
        return self._Pot(R,0.,t=t,use_physical=False)
1207919
def mass_3d(self, R, Rs, rho0, r_core):
        

        Rs = float(Rs)
        b = r_core * Rs ** -1
        c = R * Rs ** -1

        M0 = 4*np.pi*Rs**3 * rho0

        return M0 * (1+b**2) ** -1 * (0.5*np.log(1+c**2) + b**2*np.log(c*b**-1 + 1) - b*np.arctan(c))
887730
def mass_3d(self, R, Rs, rho0):
        
        Rs = float(Rs)
        m_3d = 4. * np.pi * rho0 * Rs**3 *(np.log((Rs + R)/Rs) - R/(Rs + R))
        return m_3d
859592
def dot(self, rhs):
    

    return self.x * rhs.x + self.y * rhs.y + self.z * rhs.z
102025
def ρ(self, e):
        
        d = self._dichalcogenide
        at, λ = d.at, d.λ
        return abs(2 * e - λ) * (2 * at**2)**(-1)
16790
def pot_ana(r, rho):
    
    I = 1.0
    sigma = 1.0 / rho
    phi = np.divide(I, (2.0 * np.pi * sigma * r))
    return phi
302153
def funTransSectY(theta, rho):
    
    return np.matrix([[1, rho * theta], [0, 1]], dtype=np.double)
1131953
def ystep(self):
        
        """

        self.Y[..., 0:-1] = sp.prox_l2(
            self.AX[..., 0:-1] + self.U[..., 0:-1],
            (self.lmbda/self.rho)*self.Wtvna, axis=self.saxes)
        self.Y[..., -1] = sp.prox_l1(
            self.AX[..., -1] + self.U[..., -1] - self.S,
            (1.0/self.rho)*self.Wdf)
132513
def mass_3d(self, r, rho0, gamma):
        
        mass_3d = 4 * np.pi * rho0 /(-gamma + 3) * r ** (-gamma + 3)
        return mass_3d
1085691
def _get_rho(self, w, arr):
        
        lsb = 0
        while not (w & arr[lsb]):
            lsb += 1
        return lsb + 1
520213
def _dr_dx(self, x, y):
        

        r = np.sqrt(x**2 + y**2)
        if isinstance(r, int) or isinstance(r, float):
            if r == 0:
                r = 1
        else:
             r[r == 0] = 1
        return x/r, y/r