def ncf2ffi1001(f, outpath, mode = 'w'):
    outfile = open(outpath, mode)
    header_keys = "PI_CONTACT_INFO PLATFORM LOCATION ASSOCIATED_DATA INSTRUMENT_INFO DATA_INFO UNCERTAINTY ULOD_FLAG ULOD_VALUE LLOD_FLAG LLOD_VALUE DM_CONTACT_INFO PROJECT_INFO STIPULATIONS_ON_USE OTHER_COMMENTS REVISION".split()
    print('%d, %d' % (len(f.ncattrs()) + len(f.variables), 1001), file = outfile)
    print(getattr(f, 'PI_NAME', 'Unknown'), file = outfile)
    print(getattr(f, 'ORGANIZATION_NAME', 'Unknown'), file = outfile)
    print(getattr(f, 'SOURCE_DESCRIPTION', 'Unknown'), file = outfile)
    print(getattr(f, 'VOLUME_INFO', 'Unknown'), file = outfile)
    print(f.SDATE, f.WDATE, file = outfile)
    print(f.TIME_INTERVAL, file = outfile)
    print(f.INDEPENDENT_VARIABLE, file = outfile)
    print('%d' % len(f.variables), file = outfile)
    for key, var in f.variables.items():
        print('%s, %s' % (key, getattr(var, 'units', 'unknown')), file = outfile)
    
    print(len(f.ncattrs()), file = outfile)
    for key in f.ncattrs():
        print('%s: %s' % (key, getattr(f, key, '')), file = outfile)
    
    vals = [filled(f.variables[f.INDEPENDENT_VARIABLE][:]).ravel()]
    keys = [f.INDEPENDENT_VARIABLE]
    for key, var in f.variables.items():
        if key == f.INDEPENDENT_VARIABLE: continue
        keys.append(key)
        vals.append(filled(var[:]).ravel())
        
    print(', '.join(keys), file = outfile)
    for row in array(vals).T:
        row.tofile(outfile, format = '%.6e', sep = ', ')
        print('', file = outfile)

 def compute(self, dataset_pool):
     constants = dataset_pool.get_dataset('constants')
     footprint = constants["FOOTPRINT"]
     lct = ma.filled(self.get_dataset().get_2d_attribute(self.land_cover_types), 0)
     temp = equal(lct, constants[self.lct_type.upper()])
     values = correlate(temp.astype(int32), footprint, mode="reflect")
     return self.get_dataset().flatten_by_id(ma.filled(values, 0))

 def reduce (self, target, axis=0):
     """Reduce target along the given axis."""
     axes, attributes, id, grid = _extractMetadata(target, omit=axis)
     a = _makeMaskedArg(target)
     m = getmask(a)
     if m is nomask:
         t = filled(a)
         return masked_array (numpy.maximum.reduce (t, axis))
     else:
         t = numpy.maximum.reduce(filled(a, numpy.ma.maximum_fill_value(a)), axis)
         m = numpy.logical_and.reduce(m, axis)
         return TransientVariable(t, mask=m, fill_value=fill_value(a),
                     axes = axes, grid=grid, id=id)

 def compute(self, dataset_pool):
     constants = dataset_pool.get_dataset('constants')
     fs = ma.filled(self.get_dataset().get_2d_attribute(self.footprint_size).astype(float32), 0)
     lct = ma.filled(self.get_dataset().get_2d_attribute(self.land_cover_type), 0)
     x = zeros(shape=lct.shape, dtype=float32)
     max_type = int(maximum.reduce(ravel(lct)))
     for itype in range(1, max_type+1):
         temp = equal(lct, itype).astype(int32)
         summed = correlate(ma.filled(temp, 0.0),
                            constants['FOOTPRINT'],
                            mode="reflect")
         x += temp * ma.filled(summed / ma.masked_where(fs==0, fs), 0.0)
     return self.get_dataset().flatten_by_id(arcsin(sqrt(x)))

def calcDeltaMags():
    '''
        Creates an array with the delta magnitudes of WDS objects.
    '''
    # Calculating magnitude difference for each object
    # Make nice formats...
    Pri_mag_g = ma.filled(wdsInteresting[priMag],[-999])
    Sec_mag_g = ma.filled(wdsInteresting[secMag],[-999])
    Pri_mag_gg = np.asarray(Pri_mag_g, dtype = 'float_' )
    Sec_mag_gg = np.asarray(Sec_mag_g, dtype = 'float_' )
    # Then calculate the actual delta mags
    Delta_mag = np.subtract(Pri_mag_gg, Sec_mag_gg)
    
    return Delta_mag

 def outer (self, a, b):
     "Return the function applied to the outer product of a and b."
     a = _makeMaskedArg(a)
     b = _makeMaskedArg(b)
     ma = getmask(a)
     mb = getmask(b)
     if ma is nomask and mb is nomask:
         m = None
     else:
         ma = getmaskarray(a)
         mb = getmaskarray(b)
         m = logical_or.outer(ma, mb)
     d = numpy.maximum.outer(filled(a), filled(b))
     return TransientVariable(d, mask=m)

    def create_gap_mask(self, wave):
        """
        Use the gap configuration and a wavelength vector, wave, to build a masked array that
        masks out the location of the gap.  Wavelengths between and including both gap endpoints
        will be masked out.

        Parameters
        ----------
        wave: numpy.ndarray
            Wavelength vector to construct mask from

        Returns
        -------
        mask: numpy.ma 1D masked array
            1D array masked at the locations within the configured detector gap and 1.0 elsewhere
        """
        disperser = self.instrument['disperser']
        aperture = self.instrument['aperture']
        filt = self.instrument['filter']
        gap = self.aperture_config[aperture]['gap'][disperser]
        if filt in gap:
            gap_start = gap[filt]['gap_start']
            gap_end = gap[filt]['gap_end']
        else:
            gap_start = gap['gap_start']
            gap_end = gap['gap_end']

        if gap_start is not None and gap_end is not None:
            masked_wave = ma.masked_inside(wave, gap_start, gap_end)
            mask = masked_wave / wave
            mask = ma.filled(mask, 0.0)
        else:
            mask = 1.0
        return mask

    def identify(cls, variables, ignore=None, target=None, warn=True, monotonic=False):
        result = {}
        ignore, target = cls._identify_common(variables, ignore, target)

        # Identify all CF coordinate variables.
        for nc_var_name, nc_var in target.iteritems():
            if nc_var_name in ignore:
                continue
            # String variables can't be coordinates
            if np.issubdtype(nc_var.dtype, np.str):
                continue
            # Restrict to one-dimensional with name as dimension OR zero-dimensional scalar
            if not ((nc_var.ndim == 1 and nc_var_name in nc_var.dimensions) or (nc_var.ndim == 0)):
                continue
            # Restrict to monotonic?
            if monotonic:
                data = nc_var[:]
                # Gracefully fill a masked coordinate.
                if ma.isMaskedArray(data):
                    data = ma.filled(data)
                if nc_var.shape == () or nc_var.shape == (1,) or iris.util.monotonic(data):
                    result[nc_var_name] = CFCoordinateVariable(nc_var_name, nc_var)
            else:
                result[nc_var_name] = CFCoordinateVariable(nc_var_name, nc_var)

        return result

def mag_to_flux(mag, limMag=False):
    """
    Converts magnitudes to fluxes using the following law (from Kessler+2010):
        flux = 10^(-0.4 * m + 11)
    If the magnitude is above the limit of the instrument, returns the
    corresponding limiting flux.

    INPUT:
        mag: numpy array of magnitude
        limMag: specifies the limiting magnitude

    OUTPUT:
        flux: flux-converted magnitudes
    """

    if limMag:
        exponent = (-0.4 * limMag) + 11
        limFlux = np.power(10, exponent)

        # masked arrays if magnitude is grater then limit
        maMag = ma.masked_where(mag > limMag, mag)
        maExponent = (-0.4 * maMag) + 11
        maFlux =  np.power(10, exponent)
        maFlux.set_fill_value(limFlux)

        flux = ma.filled(maFlux)
    else:
        exponent = (-0.4 * mag) + 11
        flux = 10**exponent

    return flux

 def test_testMasked(self):
     # Test of masked element
     xx = arange(6)
     xx[1] = masked
     self.assertTrue(str(masked) == '--')
     self.assertTrue(xx[1] is masked)
     self.assertEqual(filled(xx[1], 0), 0)

    def __init__(self, MetricTable, opts):

        epsilon = 0.0
        if opts.cpu:
            epsilon = 0.01

        # Create empty ratio table
        nprobs = MetricTable.nprobs
        nsolvs = MetricTable.nsolvs
        self.ratios = ma.zeros((nprobs, nsolvs), dtype=numpy.float)

        # Compute best relative performance ratios across
        # solvers for each problem
        for prob in range(nprobs):
            metrics  = MetricTable.prob_mets(prob) + epsilon
            best_met = ma.minimum(metrics)
            self.ratios[prob,:] = metrics * (1.0 / best_met)

        # Sort each solvers performance ratios
        for solv in range(nsolvs):
            self.ratios[:,solv] = ma.sort(self.ratios[:,solv])

        # Compute largest ratio and use to replace failure entries
        self.maxrat = ma.maximum(self.ratios)
        self.ratios = ma.filled(self.ratios, 10 * self.maxrat)

 def compute(self,  dataset_pool):
     buildings = self.get_dataset()  
     non_residential_sqft = buildings.get_attribute("non_residential_sqft")
     building_sqft_per_job = buildings.get_attribute("building_sqft_per_job")
     remaining_space = clip(non_residential_sqft - buildings.get_attribute("occupied_building_sqft_by_jobs"),
                            0, non_residential_sqft)
     return ma.filled(ma.masked_where(building_sqft_per_job==0, remaining_space / building_sqft_per_job), 0)

 def compute(self, dataset_pool):
     nj = self.get_dataset().get_attribute(self.number_of_industrial_jobs_wwd)
     sqft = self.get_dataset().get_attribute(self.industrial_sqft_within_walking_distance)
     regional_average = self.get_dataset().get_attribute(self.industrial_sqft).sum() / float(
         self.get_dataset().get_attribute(self.number_of_industrial_jobs).sum()
     )
     return where(sqft < 5000, regional_average, ma.filled(sqft / ma.masked_where(nj == 0, nj.astype(float32)), 0))

 def __setslice__(self, i, j, value):
     "Sets the slice described by [i,j] to `value`."
     _localdict = self.__dict__
     d = self._data
     m = _localdict['_fieldmask']
     names = self.dtype.names
     if value is masked:
         for n in names:
             m[i:j][n] = True
     elif not self._hardmask:
         fval = filled(value)
         mval = getmaskarray(value)
         for n in names:
             d[n][i:j] = fval
             m[n][i:j] = mval
     else:
         mindx = getmaskarray(self)[i:j]
         dval = np.asarray(value)
         valmask = getmask(value)
         if valmask is nomask:
             for n in names:
                 mval = mask_or(m[n][i:j], valmask)
                 d[n][i:j][~mval] = value
         elif valmask.size > 1:
             for n in names:
                 mval = mask_or(m[n][i:j], valmask)
                 d[n][i:j][~mval] = dval[~mval]
                 m[n][i:j] = mask_or(m[n][i:j], mval)
     self._fieldmask = m

    def __init__(self, MetricTable):

        # Create empty ratio table
        nprobs = MetricTable.nprobs
        nsolvs = MetricTable.nsolvs
        self.ratios = ma.masked_array(1.0 * ma.zeros((nprobs+1, nsolvs)))

        # Compute best relative performance ratios across
        # solvers for each problem
        for prob in range(nprobs):
            metrics  = MetricTable.prob_mets(prob)
            best_met = ma.minimum(metrics)
	    if (ma.count(metrics)==nsolvs and
                ma.maximum(metrics)<=opts.minlimit):
                self.ratios[prob+1,:] = 1.0;
	    else:
                self.ratios[prob+1,:] = metrics * (1.0 / best_met)

        # Sort each solvers performance ratios
        for solv in range(nsolvs):
            self.ratios[:,solv] = ma.sort(self.ratios[:,solv])

        # Compute largest ratio and use to replace failures entries
        self.maxrat = ma.maximum(self.ratios)
        self.ratios = ma.filled(self.ratios, 1.01 * self.maxrat)

 def compute(self, dataset_pool):
     hh_wwd = self.get_dataset().get_attribute(self.number_of_households_within_walking_distance)
     return 100.0 * ma.filled(
         self.get_dataset().get_attribute(self.number_of_young_households_within_walking_distance)
         / ma.masked_where(hh_wwd == 0, hh_wwd.astype(float32)),
         0.0,
     )

    def run_chunk(self, agents_index, agent_set, specification, coefficients):

        # unplaced agents in agents_index
        location_id_name = self.choice_set.get_id_name()[0]
        agent_set.set_values_of_one_attribute(location_id_name, resize(array([-1]), agents_index.size), 
                                              agents_index)
            
        ## capacity may need to be re-computed for every chunk
        if self.compute_capacity_flag:
            self.capacity = ma.filled(self.determine_capacity(capacity_string=self.run_config.get("capacity_string", None), 
                                                              agent_set=agent_set, 
                                                              agents_index=agents_index), 
                                      0.0)
            if self.capacity is not None:
                logger.log_status("Available capacity: %s units." % self.capacity.sum())
        self.run_config.merge({"capacity":self.capacity})
        if self.run_config.get("agent_units_string", None):
            self.run_config["agent_units_all"] = agent_set.get_attribute_by_index(self.run_config["agent_units_string"], agents_index)

        choices = ChoiceModel.run_chunk(self, agents_index, agent_set, specification, coefficients)

        ## this is done in choice_model
        #modify locations
        #agent_set.set_values_of_one_attribute(location_id_name, choices, agents_index)
        
        if self.run_config.has_key("capacity"):
            del self.run_config["capacity"]
            
        return choices

def flux_to_mag(flux, limFlux=False):
    """
    Converts fluxes to magnitudes using the following law (from Kessler+2010):
        flux = 10^(-0.4 * m + 11) => m = -2.5 * (log_10(flux) - 11)
    If the flux is below the limit of the instrument, returns the
    corresponding limiting magnitude.

    INPUT:
        flux: numpy array of fluxes
        limFlux: specifies the limiting flux

    OUTPUT:
        mag: magnitude-converted fluxes
    """

    if limFlux:
        limMag = -2.5 * (-11 + np.log10(limFlux))

        # applying the mask to detection below the limiting flux
        maFlux = ma.masked_where(flux < limFlux, flux)

        # to avoid warnings due to values passed to np.log10
        # fluxMask = maFlux.mask
        # maMag = -2.5 * (-11.0 + np.log10(ma.filled(maFlux,1)))
        maMag = -2.5 * (-11.0 + np.log10(maFlux))

        mag = ma.filled(maMag, limMag)
    else:
        if flux > 0:
            mag = -2.5 * (-11 + np.log10(flux))
        else:
            mag = None

    return mag

 def test_testMasked(self):
     # Test of masked element
     xx = arange(6)
     xx[1] = masked
     assert_(str(masked) == '--')
     assert_(xx[1] is masked)
     assert_equal(filled(xx[1], 0), 0)

def compute_lambda_for_vacancy(grouping_location_set, units_variable, vacant_units_variable, movers_variable,
                   secondary_residence_variable=None, multiplicator=1.0, resources=None):
        grouping_location_set.compute_variables([units_variable, vacant_units_variable, movers_variable], resources=resources)
        if secondary_residence_variable is not None:
            grouping_location_set.compute_variables([secondary_residence_variable], resources=resources)
            unitssecondary = grouping_location_set.get_attribute(secondary_residence_variable)
        else:
            unitssecondary = zeros(grouping_location_set.size())
        unitsvacant = grouping_location_set.get_attribute(vacant_units_variable)
        units = grouping_location_set.get_attribute(units_variable)
        movers = grouping_location_set.get_attribute(movers_variable)
        tsv = units - unitssecondary - unitsvacant
        lambda_value = ma.filled((units - unitssecondary).astype(float32)/ ma.masked_where(tsv==0, tsv),0) \
                            - ma.filled(unitsvacant.astype(float32) / ma.masked_where(movers==0, movers),0)
        lambda_value = lambda_value * multiplicator
        return lambda_value